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Abstract 

We propose a new method for discretizing the time variable in 
integrable lattice systems while maintaining the locality of the equa- 
tions of motion. The method is based on the zero-cm'vature (Lax pair) 
representation and the lowest-order "conservation laws". In contrast 
to the pioneering work of Ablowitz and Ladik, our method allows 
the auxiliary dependent variables appearing in the stage of time dis- 
cretization to be expressed locally in terms of the original dependent 
variables. The time-discretized lattice systems have the same set of 
conserved quantities and the same structures of the solutions as the 
continuous-time lattice systems; only the time evolution of the pa- 
rameters in the solutions that correspond to the angle variables is dis- 
cretized. The effectiveness of our method is illustrated using examples 
such as the Toda lattice, the Volterra lattice, the modified Volterra 
lattice, the Ablowitz-Ladik lattice (an integrable semi-discrete non- 
linear Schrodinger system), and the lattice Heisenberg ferromagnet 
model. For the Volterra lattice and modified Volterra lattice, we also 
present their ultradiscrete analogues. 
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1 Introduction 



The quest for a finite-difference analogue of a given differential equation can 
be justified for several sound reasons. A suitable discretization can reproduce 
most of the important properties of the differential equation in the small- 
value range of the difference interval and can be considered a "generalization" 
of the original continuous equation. Such a discretization facilitates a better 
and more intuitive understanding of the differential equation without using 
a limiting procedure, which is needed to define differentiation, and is ideal 
for performing numerical experiments. 

The suitable discretization of a completely integrable system is usually 
required to retain the integrability; if this is satisfied, it is called an integrable 
discretization. An integrable system often admits more than one integrable 
discretization; in such a case, we can consider the properties of each integrable 
discretization other than integrability and discuss which one is the most 
favorable for our purpose. The problem of integrable discretization has been 
sporadically studied since the mid-1970s, i.e., the dawn of the modern theory 
of integrable systems. For more than thirty years, various techniques have 
been developed to obtain integrable discretizations of continuous systems. 
Readers interested in the history of integrable discretizations are referred to 
the preface of Suris's book [1] . 

Partial differential equations (PDEs) involve more than one independent 
variable. The discretization of an integrable nonlinear PDE is generally per- 
formed in two steps; in the first step, we discrctize the spatial variable(s) and 
in the second step, the time variable is discretized. Of course, for some PDEs 
such as the sine-Gordon equation u^t = sin u, the roles of the space and time 
variables can be swapped and it is not meaningful to discuss the order of 
discretization. However, for an integrable nonlinear PDE wherein the roles 
of individual variables are essentially different and not interchangeable, the 
order in which the independent variables may be discretized appears to be 
unique. Thus, we can consider that the time variable is always discretized 
last, after the spatial variable(s) have been discretized. 

In this paper we focus on the problem of the integrable full discretiza- 
tion of differential-difference equations in 1 -|- 1 dimensions. The continuous 
independent variable to be discretized is regarded as time, as noted above. 
Most of the differential-difference equations considered reduce to integrable 
PDEs in a proper continuous limit; however, this is not necessary and we 
can also start with integrable lattice systems that have no continuous coun- 
terpart. The problem of time discretization is, by its nature, distinct from 
the problem of space discretization. In fact, the former problem has its own 
peculiarities and difficulties that the latter does not have. This point was 
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uncovered by Ablowitz and Ladik [2, 3] in their attempt to fully discretize 
the space-discretized nonlinear Schrodinger (semi-discrete NLS) equation [4] . 
It turned out that unexpected nonlocality emerges in the stage of time dis- 
cretization; the fully discrete NLS equation involves infinite sums and/or 
infinite products (in the case of an infinite chain) with respect to the discrete 
spatial variable and is thus a global-in-space scheme. The fully discrete NLS 
equation can superficially be written in a local form using additional depen- 
dent variables called auxiliary variables, but it does not provide any essential 
resolution of the nonlocality problem. The subsequent paper by Taha and 
Ablowitz [5] reinforces the impression that the appearance of global terms is 
a general feature of the problem of time discretization for integrable lattice 
systems. 

The pioneering work of Ablowitz and coworkers [2,3,5] is based on the 
zero-curvature (Lax pair) representation; the guiding principle is that the 
time discretization does not change the spatial part of the Lax pair for an 
integrable lattice system. This automatically guarantees the major advan- 
tages of the full discretization; that is, the time-discretized lattice systems 
have the same integrals of motion and the same structures of the solutions 
as the original continuous-time lattice systems. In more modern terms, each 
of their time discretizations belongs to the same integrable hierarchy as the 
underlying continuous-time system [6,7]. Despite the elegance of this result, 
the appearance of infinite sums/products is a shortcoming not acceptable 
to everybody, and new ideas are needed to remove the nonlocality. In this 
regard, Suris recently introduced the notion of localizing changes of vari- 
ables, applied it to a large number of integrable lattice systems, and obtained 
their time discretizations written in local equations of motion [1] . Note that 
some of his discretizations coincide with earfier results obtained using Hi- 
rota's bilinear method [8], as described in the bibliographical remarks in his 
book [1]. Suris's idea was to find, by guessing, a change of variables such that 
the equations of motion as well as the auxiliary variables can be expressed 
locally in terms of the new dependent variables; the change of variables in- 
volves the step size of time as a parameter and is considered as a discrete 
Miura transformation giving a one-parameter deformation of the original lat- 
tice hierarchy. Thus, the integrals of motion and the solution formulas for 
the time-discretized lattice system are deformed accordingly. Moreover, al- 
though Suris's approach has successfully provided many interesting examples, 
its applicability is rather limited. In particular, it is not applicable to the 
time discretization of NLS-type lattices, such as the semi-discrete NLS equa- 
tion (also called the Ablowitz-Ladik lattice) [4] , wherein the two dependent 
variables can be related by a complex conjugacy reduction. 

The main objective of this paper is to propose a systematic method for 



4 



constructing time discretizations of integrable lattice systems written as local 
equations of motion. In contrast to the other known methods, our method 
generally requires no ad hoc treatment on a case-by-case basis and appears to 
have no serious limitations in its applicability; it can be applied to possibly 
all lattice systems in 1 + 1 dimensions possessing a Lax pair representation. 
In particular, it can be used to obtain local full discretizations of NLS-type 
lattices, including the Ablowitz-Ladik lattice [4]. Actually, our method can 
be considered as a completed version of the Ablowitz-Ladik approach [2, 
3]; it both refines and extends their work in an essential way. A decisive 
breakthrough has been made by considering the lowest-order "conservation 
laws" , derived from the zero-curvature condition written in matrix form. In 
the process, a critical role is played by an arbitrary parameter in the Lax 
pair, called the spectral parameter. The requirement that all the fluxes 
corresponding to the same conserved density have to essentially coincide 
results in an "ultralocal" algebraic system for the auxiliary variables; the 
simpler case where the conserved density is trivially a constant can be treated 
in a similar manner. Thus, by solving this algebraic system, we can restore 
the locahty of the equations; that is, the global terms appearing in the stage 
of time discretization can be replaced by local expressions in terms of the 
original dependent variables. 

This paper is organized as follows. In section 2, we describe the general 
method for discretizing the time variable. In section 3, we construct time 
discretizations of the Toda lattice in Flaschka-Manakov coordinates [9-11], 
the Ablowitz-Ladik lattice, the Volterra lattice, the modifled Volterra lat- 
tice, and the lattice Heisenberg ferromagnet model. As a spin-off, we obtain 
ultradiscrete analogues [12, 13] of the Volterra lattice and modified Volterra 
lattice. In addition, we uncover unexpected relationships with the work of 
Nijhoff, Quispel, Capel et al. [14-17]. Section 4 is devoted to concluding 
remarks. 

2 Method for time discretization 

In this section, we discuss the problem of time discretization for a given inte- 
grable lattice system. We start with a Lax pair formulation in the continuous- 
time case and then proceed to discretize the time variable. A set of auxiliary 
variables is introduced to express the time-discretized lattice as a closed sys- 
tem of equations. Using the fundamental "conservation laws" derived from 
the Lax pair, we can obtain local expressions for the auxiliary variables in 
terms of the original variables. 
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2.1 Lax pair and a conservation law in the continuous- 
time case 

The Lax pair formulation in a semi-discrete space-time comprises a pair of 
linear equations, 

= Ln{X)^n, ^n,t = M„(A)^„. (2.1) 

Here, \E'„ is a column-vector function and n is the discrete spatial variable. 
The subscript t denotes differentiation with respect to the continuous time 
variable t. The square matrices L„ and M„ depend on the spectral parameter 
A, which is an arbitrary constant independent of n and t. The compatibility 
condition of the overdetermined system (12.11) is given by [2,4,6,7,18,19] 

Ln,t + LnMn - M„+iL„ = O, (2.2) 

which is (a semi-discrete version of) the zero- curvature condition. The sym- 
bol O is used to stress that this is a matrix equation. If we specify the 
A-dependent matrices L„ and M„ appropriately, (12.21) results in a closed 
differential-difference system for some A-independent quantities in L„ and 
Mn. In such a case, the pair of matrices L„ and M„ is called a Lax pair. 
The matrix L„ is usually ultralocal in the dependent variables; that is, if 
Ln involves some variable, say then it does not involve shifted variables 
such as Un±i and Un±2- In addition, the determinant of L„ is required to be 
nonzero for generic A so that the spectral problem is well-posed on the entire 
infinite chain. 

The zero-curvature condition (12.21) generates a conservation law of the 
following form: 

^log(detL„) = A+(trM„). (2.3) 
Here, is the forward difference operator in the spatial direction, i.e., 

Kfn := Ui - fn. (2.4) 

For a proper integrable lattice system, detL„(A) has to be either a time- 
independent function of A or the exponential of a A-independent conserved 
density multiplied by an overall t-independent factor. Indeed, if this was not 
satisfied, e.g., detL„(A) = 1 -|- Xun, then the expansion of log det L„(A) with 
respect to A would yield an infinite number of (almost) ultralocal conserved 
densities, say m„, u^, . . . , and thus the lattice system would be trivial in 
some sense. 

Note that the zero-curvature condition (12. 2p has the following invariance 
properties: 
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(a) Ln /(A)L„, where /(A) is a t-independent scalar function, 

(b) M„ M„ + g{X)I, where g{X) is an n-independent scalar function and 
/ is the identity matrix, 

(c) Ln — * exp{apn)Ln and M„ —>■ Mn + ajnl, where a is a parameter, p„ 
is a conserved density and jn is the corresponding flux (up to a sign), 
namely, dtPn = A+j„. 

In particular, using f l2.3p and the above properties, we can convert the Lax 
pair to a normalized form, i.e., detL„ = 1 and trM„ = 0. 

2.2 Lax pair in the discrete- time case and auxiliary 
variables 

Now, we discuss how to construct the time discretization of a given lattice 
system having the Lax pair Ln and Mn- The natural discrete-time analogue 
of the linear system (12.11) is given by 

^„+l = LniX)^n, $n = K(A)^n, (2.5) 

where the tilde denotes the forward shift (m — * m + 1) in the discrete time 
coordinate m G Z. Here and hereafter, the dependence on m is usually sup- 
pressed unless it is shifted. The compatibility condition of this overdeter- 
mined linear system is given by [2,3,5-7,19,20] 

LnVn = Vn+lLn, (2.6) 

which is (a fully discretized version of) the zero- curvature condition. Note 
that (12. 6p can be rewritten as 

Vn = L-^Vn+lLn OT LnV^L-' = K+l- (2.7) 

Following the work of Ablowitz and coworkers [2,3,5], the matrix L„(A) is 
assumed to be the same as that in the semi-discrete case. Then, we look 
for a Vn{X) such that the zero- curvature condition (12. 6p results in a closed 
system of partial difference equations providing a discrete-time analogue of 
the semi-discrete system. 

For this purpose, we assume that the matrix Vn has asymptotic behavior, 

Vn = I + h[Mn + Oih)], (2.8) 

where h is a sufficiently small (but nonzero) parameter independent of A and 
is considered the difference interval of time (cf. (12. 2p and (12. 6p ). More pre- 
cisely, h approximates the "true" step size of discrete time up to an o{h) error. 



7 



Moreover, we assume that Ki(A) has essentially the same A-dependence as 
J + hMn{X). Note, however, that M„(A) can only be determined up to the 
addition of an n-independent scalar matrix (cf. (b)). Thus, this arbitrariness 
has to be taken into account in determining Vn{\); this corresponds to the 
freedom (nonuniqueness) of choosing the linear part of the time-discretized 
lattice system that determines the dispersion relation [3]. Moreover, some 
n-independent quantities (usually set as constants) in M„ translate into n- 
dependent quantities in Vn, which typically constitute, up to the reformula- 
tion of the dependent variables, new auxiliary variables. Then, we specify 
appropriate boundary conditions for these new variables in Vn, which should 
retrieve the corresponding n-independent values in M„. In fact, we usually 
assume "constant" boundary conditions for Vn, 

lim Vn = lim Vn = n-independent and finite. (2-9) 

n— oo n— >+oo 

However, the right-hand side is allowed to depend on the time variable m G Z. 
In the application of the inverse scattering method based on the Lax pair, 
we need to specify the boundary conditions for L„ as — > ±00. In such a 
case, it is redundant to impose the boundary conditions on Vn at both spatial 
ends as given in fl2.9p . and it is nontrivial that the redundant boundary 
conditions are compatible. In fact, it is sufficient to know only one of the two 
boundary values, lim„^_oo Vn or lim„^+oo Vn- In the existing literature [1-3, 
6,7,19], it is hypothesized that these two limits indeed coincide; a preliminary 
consideration without using this hypothesis is given in [21]. In section [HI we 
demonstrate for specific examples that this is not a hypothesis but a verifiable 
fact. 

Let us decompose L„(A) into a sum of terms, each of which is the product 
of an (n, m) -independent scalar function of A and a A-independent matrix, 
i.e., 

^max 

^n(A) = f^^^^^n- 

1 — ^min 

The scalar functions /i(A) (imin < < ^max) are linearly independent; typi- 
cally, they are powers of A, e.g., fi{X) = A*. The nonzero elements of the 
matrices Ln^ are classified into two types, that is, constants (or, at most, 
arbitrary functions of only one independent variable) and dynamical vari- 
ables depending on both independent variables. We express the entire set 
of functionally independent dynamical variables in Ln^ {imin < < ^max) as 
{l„}. The set of dynamical variables {vn} is defined from Ki(A) in exactly 
the same way. Then, the zero- curvature condition (12. 6p provides a (typically 
bilinear algebraic) system for {Z„, Z„} and {Vn, Vn+i}- In particular, this sys- 
tem contains a useful subsystem, that is, an ultralocal and linear system in 
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{Vn}, wherein the coefficients involve and its shifts; this subsystem can 
also be derived by noting that fl2.7p holds true as an identity in A. We solve 
this subsystem to express a subset of {vn} in terms of the remaining {vn} 
as well as {/„} and its shifts such as {/„} and {/„-i}. Thus, we can reduce 
the number of independent dynamical variables in Vn while maintaining the 
ultralocality of Vn with respect to {vn}- It remains to be verified that the 
zero-curvature condition (12.61) indeed provides a meaningful fully discrete 
system for the reduced set of dependent variables; it should define a con- 
sistent and unique time evolution for generic initial data under appropriate 
boundary conditions (cf. fl2.9l) ). 

2.3 Fundamental "conservation laws" in the discrete- 
time case 

We consider the determinant of both sides of (12. 6p to obtain the equality 

(detL„)(defl4) = (det L„)(det K+i)- (2.10) 

This relation can be written more explicitly in the form of a discrete conser- 
vation law, 

A+log(detL„) = A+log(detK). 

Here, is the forward difference operator in the time direction, i.e., A^/„ := /, 
This is the discrete-time version of (12.31) . Note that after cancelling the m- 
independent factor of det L„, the relation (12.101) should reduce to either 

det K, = det K+i (2.11) 

or 

exp [pn) det Vn = exp (p„) det K+i, (2.12) 

where pn is a nontrivial conserved density. In the ffist case (12.111) . one may 
speculate that the above conservation law may become the meaningless re- 
lation A^ (const.) = A^ (const.). However, in all the examples that require 
the introduction of auxiliary variables, this appears not to be the case; it is 
not immediately evident that the determinant of Vn is n- independent, and 
thus the relation (I2.10p still contains meaningful information. We employ a 
simplified but still ultralocal (with respect to {f «}) form of Vn, compute its 
determinant, and expand it with respect to A in the summed form 

J max 

detV;(A)= Yl 9jW^- (2.13) 

J— Jmin 
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Here, gj{X) (jmin < J < linearly independent functions of A, e.g., 

gj{X) = X'^\ and their coefficients a^^ are A-independent functions of 
{In}, and the space/time shifts of {Z„}. Substituting f l2.13p into (12. lip , we 
obtain the n-independence of the coefficients of gj{X) for all j, i.e., 



a. 



(f)= hm ali\ Jmin<J<Jmax. (2.14) 



^±00 



In the second case f l2.12p . which is more common than (12.111) . the substitution 
of (I2.13P gives the set of relations 

exp (pn) a^f ^ = exp (p„) ^^j^, j^m <j< jmax- 

Thus, there exist seemingly more than one flux logOn '* associated with the 
same conserved density p„, but the an^ should coincide up to trivial pro- 
portionality factors. Indeed, calculating the ratio of the above equality for 
different values of j on both sides, we obtain 

~Th) = ^^"^ ~Th)^ ^^^^ - ^1 ^ ^2 < Jmax- (2.15) 

Lin iJ'n 

That is, the ratio ai^"""'' : ■ ■ ■ : an™'''''' is independent of n. In both the above 
cases, the right-hand side of (I2.14p or (I2.15P is determined by the bound- 
ary conditions for Vn, in particular, the boundary values of the dependent 
variables contained in Vn] each right-hand side is set as a definite value in- 
dependent of n (cf. (12. 9p ). This results in an "ultralocal" algebraic system 
for a subset of {vn} that essentially constitutes the auxiliary variables; the 
number of independent unknowns is usually equal to that of the indepen- 
dent equalities so that this algebraic system is neither overdetermined nor 
underdetermined. 



2.4 Algebraic system for the auxiliary variables and 
local equations of motion 

We solve the obtained algebraic system for the auxiliary variables appearing 
in Vn (and not in M„). By eliminating all but one of the auxiliary variables, 
this system becomes a scalar algebraic equation in the remaining auxiliary 
variable; the other auxiliary variables are expressible in terms of the solution 
of this equation. The degree and complexity of this algebraic equation depend 
on the boundary conditions for Vn, which determine the dispersion relation 
of the time-discretized lattice system. In this paper, we mainly consider the 
case for a degree of two so that the equation is solvable by the quadratic 
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formula. For sufficiently small h, referring to the prescribed behavior fl2.8p of 
Vn, we can discard one of the two solutions as improper, and obtain a unique 
proper solution of the quadratic equation. In general, the larger the number 
of grid points defining the lattice system as well as the dispersion relation, 
the higher the degree of the algebraic equation determining the auxiliary 
variables. The higher-degree case can be interpreted as a composition of 
lower-degree cases; that is, the time evolution in the higher-degree case can 
be factorized into sequential applications of more elementary time evolutions. 
This point will be illustrated in section [3l Once all the auxiliary variables 
have been expressed locally in terms of the original dependent variables that 
have already appeared in the semi-discrete case, we only have to substitute 
them into an appropriate subset of equations arising from the zero-curvature 
condition fl2.6p . Because of the use of the "conservation laws" (cf. (12. lip 
or (12.121) ). not all of the equations arising from (12. 6p are independent and 
necessary any longer. We choose a minimal subset of these equations so that 
the substitution of the local expressions in the auxiliary variables produces 
the unique discrete-time evolution of the lattice system. 

Last but not least, the above set of conservation laws used to determine 
the auxiliary variables can, in principle, be derived from the original system of 
partial difference equations resulting from the zero-curvature condition (12.61) . 
However, in practice, this is an extremely difficult task; any computation 
conducted at the component level without following a set procedure is highly 
unlikely to arrive at the nontrivial conservation laws to be derived. Thus, 
our derivation performed at the matrix level using a determinantal formula 
and the spectral parameter is possibly the only way of obtaining them. 

2.5 Remarks on nonautonomous extensions 

Actually, the "step size" parameter h introduced in (12. 8p need not be a con- 
stant and can depend arbitrarily on the discrete time coordinate m G Z; in 
other words, the discrete-time flow involving an arbitrary parameter h be- 
longs to the same hierarchy, and we can specify any value of h at every step 
of the time evolution. Indeed, this does not result in any essential difference 
in the subsequent computations, because the zero-curvature condition (12. 6p 
involves equal-time Vn only. In this way, we can obtain nonautonomous ex- 
tensions of time-discretized lattice systems involving one arbitrary function 
of the discrete time (see, for example, [22] and references therein). More- 
over, if the time discretization considered can be factorized into a composi- 
tion of M elementary time evolutions (cf. [6, 7, 19]), then it generally allows 
an extension involving M arbitrary functions of time; however, we do not 
proactively discuss this possibility to avoid unnecessary confusion. From the 
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modern point of view (cf. [20,23-25]), a discrete-time flow can be identi- 
fied witfi tfie spatial part of an auto-Backlund transformation of the original 
continuous-time flow or an auto-Backlund transformation of the whole hier- 
archy of continuous-time flows. Thus, the consistency of this nonautonomous 
( "nonuniform in time" ) extension can be understood as the commutativity of 
auto-Backlund transformations for different values of the Backlund param- 
eter [26,27], known as Bianchi's permutabihty theorem. Readers interested 
in Backlund transformations are referred to the Proceedings [28] . 

Bianchi's permutabihty often implies an ultralocal relation among the 
four solutions of the same system, called a nonlinear superposition formula; 
its derivation is based on the compatibility condition of two "different" 
time evolutions: ^ = V{X, hi)^ and ^ = F(A, /i2)^- If there exist two or 
more distinct one-parameter auto-Backlund transformations, then the ma- 
trix V{\,h2) may be generalized to W{\,h2). The nonlinear superposition 
formula motivates us to place the four solutions appropriately at the four ver- 
tices of a rectangle [15, 16] and to assign each value of the Backlund parameter 
to each pair of parallel sides; this is often referred to as the Lamb diagram 
(see Fig. 3 in [26]). It is convenient to identify each value of the Backlund 
parameter with the length of each side, although values are not restricted to 
positive numbers. We can use this rectangle as an elementary cell defining 
two directions of new independent variables [15-17]. In fact, repeated ap- 
plications of the auto-Backlund transformation at (generally) distinct values 
of the Backlund parameter generate a two-dimensional unequally spaced (in 
both directions) lattice in the quadrant (cf. Fig. 12 in [27] and the main figure 
in [29]), or even in the plane, with elementary cells of various sizes. With this 
understanding, we can reinterpret every nonlinear superposition formula as 
a fully discretized lattice system involving one arbitrary function of one in- 
dependent variable and another arbitrary function of the other independent 
variable. 

Thus, if we encounter a time discretization of an integrable lattice system 
that has the same form as a nonlinear superposition formula for some con- 
tinuous/discrete integrable system, then it allows a natural nonautonomous 
extension involving two arbitrary functions originating from two values of 
the Backlund parameter. For example, the nonlinear superposition formula 
for the potential KdV hierarchy (cf. (14) or (16) in [30]) suggests a nonau- 
tonomous extension of a discrete potential KdV equation (cf. (5) and (6) 
in [31]; (5.16) in [32]; (68) in [33]), while the nonlinear superposition formula 
for the sine-Gordon equation (or the potential mKdV hierarchy) with sim- 
ple sign handling implies a nonautonomous extension of the fully discretized 
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sine-Gordon equation [20, 34] given by 



ri+l,m ,m+l 1 /(-> 1 

tan ■ — — = —— ——tan ■ — ■ , (2.16) 

V 4 y f{n)-g{m) \ 4 / v / 



or equivalently, 



sm ' 



g(m) V 4 
= sm . (2.17) 



If we rescale/redefine the arbitrary functions and variables as f{n) = {4/A)F{nA) 
and g{m) = hG{mh), and Un,m = u{nA,mh), nA =: x, and mh =: t, respec- 
tively, then the continuous limit A,h ^ reduces f l2.17p to the variable- 
coefficient sine-Gordon equation 

F{x) 

-Q^iJ'xt = smw, 

which is obviously equivalent to the constant-coefficient sine-Gordon equation 
uxT = sinu. Thus, this type of nonautonomous extension in the discrete case 
is thought to be a vestige of coordinate transformations in the continuous 
case that do not mix the two independent variables. 

It is still unclear whether any of the time-discretized lattice systems ob- 
tained in the next section in their present form can be identified with a 
nonlinear superposition formula. Alternatively, we propose an intriguing 
procedure for constructing fully discrete nonautonomous systems involving 
two arbitrary functions; this procedure is closely related to the property of 
three-dimensional consistency [35-37]. It is assumed that L„(A) is ultralocal 
in the set of dynamical variables {In}, while V^(A, h) can be written in terms 
of {ln-i,ln} or {ln-i,ln}- Wc ouly cousidcr the former case, because the 
latter case can be dealt with in a similar manner. Then, the linear problem 
(12.51) with the "step size" h set as h = hi and h = h2 respectively implies the 
relations 

^„ = iy«^„_i, iV« := K(A, /ii)L„_i(A), (2.18a) 

and 

$„ = iVf iVf := K(A, /i2)^n-i(A). (2.18b) 

Here, the newly introduced matrix Nn^ connects the values of \I' at two 
lattice points (n — 1, m) and (n, m + 1), and depends only on the dynamical 

^ (2) 

variables at these points, namely, {Z„_i,Z„,}. Similarly, connects the 
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two points (n — 1, m) and (n, m + V), and depends on the dynamical vari- 
ables at these points: {In-i, In}- Note that {n, m + 1') can be identified with 
(n, m + 1) only if hi = The compatibility condition of (12.181) is given by 
a "new" zero-curvature equation, 

Thus, the substitution of the explicit forms of the matrices Nn\{ln-.i, In}, hi) 
and Nn {{ln~i, In}, ^2) into this equation should provide a closed system for 
the dynamical variables at the four lattice points 

In—l, In, In, ^n+l( In+l)- 

We can regard the corresponding parallelogram as defining two directions of 
new independent variables, so that the system now involves two arbitrary 
functions originating from hi and h2. 



3 Examples 

In this section, we apply the general method for time discretization in section 
2 to five important examples: the Toda lattice, the Ablowitz-Ladik lattice 
(semi-discrete NLS), the Volterra lattice, the modified Volterra lattice, and 
the lattice Heisenberg ferromagnet model. 

We are only interested in time discretizations that can accurately approx- 
imate the continuous-time dynamics in the small-value range of the "step 
size" parameter h. This implies that the discrete-time Lax matrix Vn al- 
lows the asymptotic expansion with respect to h given in (12. 8p : in particular, 
(yn — I)/h does not involve 0{l/h) terms. This requirement plays a cru- 
cial role in obtaining suitable local expressions for the auxiliary variables 
appearing in Vn- 

The parameter h is generally assumed to be nonzero; alternatively, one 
can allow the case h = as the trivial identity mapping Z„ = Z„. 

3.1 The Toda lattice in Flaschka— Manakov coordinates 

We consider the Toda lattice written in Flaschka-Manakov coordinates [9- 
11]: 

Un,t = Un{Vn - Vn-l), t^n,t = - M„. (3.1) 

The parametrization 

Un = e " " ^, Vn = Xn,t, (3-2) 
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enables the system (13. 11) to be rewritten as the Newtonian equations of motion 
for the Toda lattice, 



■^n,tt — " " 



The Lax pair for the Toda lattice (13. ip in Flaschka-Manakov coordinates is 
given by [1] 



-1 



n ""n 



Un 





-Un 

1 A + Vn^l 



(3.3a) 



(3.3b) 



Indeed, the substitution of (13. 3p into the zero- curvature condition (12. 2p re- 
sults in (13. ip . 

Let us move to the discrete-time case. A comprehensive overview of the 
existing results on time discretizations of the Toda lattice is given in §3.22 
and §5.11 of [1], thus we do not repeat it here. In view of the zero-curvature 
condition (12. 6p . we assume the Lax matrix Vn of the following form: 



Vn = I + h 



-Xa + On —aUn — Unbn 

a + hn Xbn + a-a-l + aVn-1 + bnVn-1 



(3.4) 



Here, a„ and 6„ are auxiliary variables. The zero-curvature condition (12. 6p 
for the Lax pair (I3.3ap and (13. 4p amounts to the following system of partial 
difference equations: 



f 1 

h 
1 

h 

an - On+l = a{Vn- Vn) , 



{Un - Un) = Otn+lUn - + Vn {aUn + Unbn) " {aUn + Un&n) Vn-U 

{Vn - Vn) = an+lVn - VnQn + a {Un+1 - Un) + Mn+l&n+l " ^n&n, 



(3.5) 

The general form (13. 5p of the time-discretized Toda lattice is integrable for 
matrix- valued dependent variables, but in the following, we consider only the 
case of scalar dependent variables. In view of (12. 9p and (13. 2p , we impose the 
following boundary conditions for m„, t>„, a„, and 6„: 



lim M„ = 1, lim Vn = 0, lim a„ = 0, <y„ 



n— »±oo 



n— >itoo 



n— »±oo 



lim bn = b. 

n— >±oo 



(3.6) 



The boundary value of a„ is set as zero by redefining h. To be precise, 
the boundary conditions (13. 6p contain redundant information. Indeed, it can 
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be shown that the auxihary variables have the same hmit values for n ^ — oo 
and n +00. Therefore, it is sufficient to assume either lim„^_oo(an5 = (0, 
or lim„_>4.oo(ttn, bn) = (0, b). In the continuum limit of time /i — 0, the aux- 
iliary variables reduce to constants, i.e., a„ ^ and 6„ b. Thus, in this 
limit, the time discretization (13.51) reduces to 

Un,t = {a + b)Un{Vn - fn-l), Vn,t = {a + b){Un+l - Mn), 

which is, if a + 6 7^ 0, equivalent to the Toda lattice (13. ip up to a redefinition 
of t. Note that in the case a + b = 0, (13. 5p has the trivial solution m„ = 

Vn = fn, "n = 0, 6„ = 6. 

The determinant of the 2x2 Lax matrix L„ (I3.3al) can be immediately 
computed as detL„ = Un- Using the recurrence formula for in (13.51) . we 
can rewrite the Lax matrix Vn (13. 4p in an ultralocal form with respect to the 
auxiliary variables and Thus, its 2 x 2 determinant is computed as 

detK(A) 

= [l + h (-Aa + an)] [I + h (A6„ + a„ + at;„_i + b^Vn^i)] 

+ {aun + uj)n) {a + 6„) 
= —X^h^abn — \h [(a — 6„) (1 + /ia„) + ha?Vn-i + habnVn-i\ 

+ (1 + ha^f + (1 + han) {havn-i + hbnVn-i) + h'^ {aUn + Unbn) {a + bn) ■ 

Therefore, the equality (12.101) combined with the boundary conditions (13. 6p 
(or the streamlined version as stated above) implies the set of relations 

( bn = bAn, 



a - 6„) (1 + han) + ha Vn~i + habnVn-i = {a - b) An, 



1 + hanf + (1 + han) {havn-i + hbnVn-l) 



(3.7) 



+ h'^ {aUn + Unbn) + &n) = [l + /i^ (o + 6)^] A„. 

Here, A„ is the quantity that satisfies the recurrence formula 

UnAn = M„A„+i, (3.8) 

and has the normalized boundary value, lim„^-|-ooAn = 1. More precisely, 
only one of the two conditions lim„^„oo A„ = 1 and lim„_»+oo A„ = 1 is re- 
quired, in accordance with the boundary conditions for a„ and 6„. Then, the 
other condition can be confirmed. Thus, A„ can be written explicitly (but 
globally) as 

n— 1 ~ +00 

Ui -r-r Ui 



An= n ^=n- 

J. J. 1/ . J. J. 7/ . 



Un ^ ^ U4 

-OO J=n 
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Note that the first equahty in fl3.7p is vahd even if a = (cf. fl3.5p ). The al- 
gebraic system fl3.7p essentially comprises two nontrivial relations containing 
the two unknowns, and A^. 

It is rather cumbersome to identify the dispersion relation of (13.51) from 
its linear leading part (cf. subsection 13.21) . Instead, we can compute the 
asymptotic form of the Lax matrix Vn (13.41) as n ^ — oo or n — >■ +00 and 
consider its factorization (cf. [1,6,7,19]). The results imply that the cases 
a = and 6 = form the basis for the general case; that is, the time evolution 
for general a and b is equivalent to the order-independent composition of two 
time evolutions corresponding to a = and 6 = 0, respectively. In addition, 
these two cases are related to each other through the time refiection, as we 
will see below. 

Now, we consider the two fundamental cases: a = or 6 = 0. 

• The case a = 0, 6 7^ 0. 

In this case, the auxiliary variable a„ vanishes. Thus, the discrete-time 
system (13. 5p reduces to (cf. (3.6) in [38]) 

' 1 _ 

- (m„ - M„) = 6A„+lM„ {Vn -Vn-l) , 

UnAn = UnAn+1, lim A„ = 1 or lim A„ = 1. 

^ n— >— oo n— >+oo 

The algebraic system (13.71) simplifies to a quadratic equation in A„, 

{hbfunAl - [1 + {hbf - hbvn^i] A„ + 1 = 0. (3.10) 

The asymptotic behavior (12.81) of the Lax matrix Vn implies that the proper 
solution of (I3.10p is given by 

[1 + [hbf] A„ = ^ ^ (3.11) 

1 - evn-i + y(l - evn-if - 4e2u„ 

where e := hb/ [1 + [hbf]. Thus, if hb G M, then -1/2 < e < 1/2. In this 
case, the local expression (13.111) is valid only if — 1 < /i6 < 1, which covers 
the range — 1/2 < e < 1/2. The borderline cases hb = ±1 are excluded from 
our consideration. If (hb)"^ > 1, (13.110 is inconsistent with the boundary con- 
ditions, and the other solution of (13.101) should be adopted. In any case, the 
boundary conditions for Un and Vn imply that lim„^_oo A„ = lim„^_|_oo A„ = 1. 
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Substituting fl3.11l) into the first and second equations in (I3.9p . we obtain 
a time discretization of tlie Toda lattice fl3.1l) in tfie local form, 

2Un (Vn -Vn-l) 



2 iUn+1 - Un) 



(3.12) 



1- evn + \ (1 - evn 



4e2u 



■n+l 



This is an explicit scheme under the boundary conditions lim„^±oo = 1 
and lim^^-too t^n = 0; Vn is determined from «„, t>„, and Un+i, and subse- 
quently, Un is determined from f„_i, m„, f„, Vn, and m„+i. If /i& G R, the 
Un are nonzero and real-valued, and the f„ are real-valued at the initial 
time, then (13.91) implies that the auxiliary variable A„ is always real-valued. 
Consequently, the discriminant of the quadratic equation (13.101) must be non- 
negative. Thus, as long as the amplitudes of m„ — 1 and f„ are sufficiently 
small and their effects can be regarded as perturbations, the real-valuedness 
of Un and Vn is preserved under the time evolution of the discrete-time Toda 
lattice flHT^ for -1/2 < e < 1/2. 

To express the backward time evolution explicitly, we only have to re- 
place A„+i in the first and second equations of (13.91) with AnUn/ Un and then 
substitute the local expression (13. lip . The resulting system is 

2Un (Vn - Vn-l) 



Ur. 



Ur 



1 

Ur, 
U„ 



eVn-1 + 



1 - evn-if - 4:e^Ur. 

2 {Un+l — Un) 



(3.13) 



1 - evn-i + \J {I - evn-i) -^e^Un 
Using the first equation, the second equation in (13.130 can also be written as 

2 (Un+l - Un) 



Vn - Vr. 



{2Vn - Vn-l) + a/(1 - (^Vn-lY 



4:6^ Un 



• The case 6 = 0, a 7^ 0. 

In this case, the auxiliary variable bn vanishes. Thus, the discrete-time 
system (13. 5p reduces to 



— (1 + han + haVn 
h 

— (1 + han + haVn 
h 



a iv„ 



^1) [Un 

(Vn - 
- V„) , 



Ur 



aUn {Vn 



'Jn-l) , 



a [Un+l 



Ur 



(3.14) 



lim Or, = or lim 



0. 
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The algebraic system fl3.7l) simplifies to 



1 + han + havn-i = A„, 

(1 + han) + (1 + han) haVn~i + (ha) Un = [1 + iha) \ A^. 

Note that, using A^i instead of system f l3.14p with (13.81) and (13.151) can 
be identified with the discrete-time Toda lattice given in §3.8 of [1] (also 
see [39]). 

By eliminating A„, (13.151) reduces to a quadratic equation in 1 + han-, 
(1 + hanf — [l + (haf — hcuVn-\\ (1 + /ia;„) + (/ia)^M„— [l + (/la)^] haVn-i = 0. 
The proper solution of this quadratic equation is given by 

1 + han ^ " '^^"-1 + \/(l + SVn-l)"^ - A5'^Un 



1 + {ha) 



(3.16) 



where 6 := ha/ [1 + {ha)"^]. Thus, if ha G M, then -1/2 < 5 < 1/2. In this 
case, the local expression (13.161) is valid only if — 1 < /la < 1, which cov- 
ers the range — 1/2 < 5 < 1/2. The borderline cases ha = ±1 are excluded 
from our consideration. If (/la)^ > 1, (I3.16P is inconsistent with the bound- 
ary conditions, and the other solution of the quadratic equation should be 
employed. In any case, the boundary conditions for u„ and f„ imply that 
lim„^__oo an = lim„^+oo = 0, and consequently, lim„^_oo A„ = lim„^+oo A„ 

Substituting the local expression (I3.16P into the first and second equations 
of (I3.14p . we obtain a time discretization of the Toda lattice (13. ip . 

Un-Un 2Un (fn " Vn-l) 



1 + SVn^^ + J(l + SVn^^)^ - ^5^Un 

\( ~^ ^3.17) 

Vn-Vn 2 [Un+l - Un) 



l + 6{2Vn-Vn-l) + \l{l + 6Vn-lY -Ad^Un 

Under the boundary conditions lim„^-|-oo = 1 and lim„^-|-oo = 0, this is 
an explicit scheme in the forward time direction; Un is determined from Wn-i, 
M„, and fn, and subsequently, Vn is determined from Wn-i, m„, m„, and 
Un+i- As long as m„ and Vn can be regarded as perturbations around the 
boundary values 1 and 0, their real-valuedness is preserved under the time 
evolution of the discrete-time Toda lattice (13.170 for —1/2 < 5 < 1/2. 

To obtain the backward time evolution explicitly, we rewrite a„ in the 
first and second equations of (I3.14p as an+i + a (vn — f„) and then substitute 
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the local expression fl3.16p . The resulting system is 



Ur 



Ur. 



2Un {Vn-Vn-l) 



2 [Un+l - Ur. 



n+1 



(3.18) 



It is easy to see that fl3.17p is equivalent to fl3.13l) through the time reflection 
and the identification 5 ^ —e. In the same manner, (13.181) can be identi- 
fied with (I3.12P through the time refiection. Thus, the forward/backward 
time evolution in the case 6 = corresponds to the backward/forward time 
evolution in the case a = 0, up to a redefinition of the parameters. 



3.2 The Ablowitz-Ladik lattice 

In this subsection, we discuss the time discretization of the Ablowitz-Ladik 
lattice, which we consider to be the most instructive example to illustrate 
our general method. The (nonreduced form of the) Ablowitz-Ladik lattice is 



qn,t - aqn+i + hqn^i + (a - h)qn + aqn+ir^qn - bqnTnqn-i = 0, 
rn,t - hvn+i + ar„_i + (6 - a)r„ + 6r„+ig„r„ - ar^qurn-i = 0, 



(3.19) 



and its Lax pair (cf. (12. ip and (12. 2p ) is given by 



A 



+ 



qn 
r„ 



+ 







A qn 
r i 



(3.20a) 



A^a - a(l + qnrn-i] 



XaVr 



Xaqn + -^qn-i 
^ - 6(1 + r„g„_i) 



(3.20b) 



Here, the free parameters a and b are usually set as constants, but they can 
depend on the time variable t in an arbitrary manner, namely, a := a{t), 
b := b{t). The familiar form of the Ablowitz-Ladik lattice is obtained from 
(I3.19P through the reduction of the complex conjugate: b = a*, r„ = ag*, 
where a is a real constant. The case where a and b are purely imaginary 
leads to an integrable semi-discretization of the NLS equation. Taking into 
account the A-dependence of the Lax matrix M„ fl3.20bp . we look for a Lax 
matrix Vn with the following dependence on A: 



VnW = i + h{ AVp) + Av;« + + \vy^ + -1^!-) 



(3.21) 
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Here, Vn^ {j = 2, 0, —2) are diagonal matrices, while Vn' {j = 1, —1) are 
off-diagonal matrices. The nonzero parameter h may depend on the discrete 
time coordinate m G Z. Substituting fl3.20ap and fl3.2ip into (12. 6p . we find 
that the matrix VC, should assume the form 



Vn = I + h 





X{aqn - qnCn) 

+ x(6g„-i - dnqn-i) 


+ \{hrn -rndn) 


y?Cn + rndnqn-1 



(3.22) 



Here, a„, c„, and dn are auxiliary variables. The zero-curvature condition 
(12.60 is equivalent to the following system of partial difference equations: 



T {qn - qn) - aqn+l + 6gn~l " "n+ign + ^n^n 

+ g„+ic„+i(l - r„g„) - (1 - qnrn)dnqn-i = 0, 



6r, 



n.+l 



+ ar, 



n-l 



+ r„+irf„+i(l - g„r„) - (1 - rnqn)cnr.. 



n—l 



(3.23) 



- ttn = a(g„r„_i - g^+i^n), 

/3n+l — Pn = b{rnqn-l — ^n+iqu), 

(1 - r„g„)c„ = c„+i(l - r„g„), 
(1 - qnrn)dn = dn+iil - g„r„). 

Actually, the general nonreduced form (13.230 of the time-discretized Ablowitz- 
Ladik lattice is integrable for matrix-valued dependent variables. However, 
in this paper, we consider only the case of scalar dependent variables. In ad- 
dition to (12. 9p . we impose rapidly decaying boundary conditions for g„ and 
r„, that is, 

lim g„ = lim r„ = 0, lim (a„, c„, (i„) = (a, /3, c, d) . (3.24) 

n— >itoo n^±oo n— >itoo 

Similarly to the semi-discrete case, the parameters a, b, a, P, c, and d are 
allowed to depend on the discrete time coordinate m G Z,. Thus, the time 
dependence of h can be absorbed by a, b, an, Pn, Cn, and dn- For the same 
reason, it is possible to set h as unity, but we prefer to leave it as a small pa- 
rameter, usually < ^ 1. To be precise, the boundary conditions (I3.24p 
contain redundant information. Indeed, it can be shown that the auxiliary 
variables have the same limit values for n —oo and n ^ +oo. There- 
fore, it is sufficient to assume either lim„^_oo {c(n, l3n,Cn,dn) = {a,P,c,d) or 
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lim„^+oo {an, Pn, Cn, d^) = (a, P, c, d). The last four relations in (I3.23P imply 
the following global expressions for the auxiliary variables c„, and dn 

in terms of g„ and r„ [1,3,5,7,19]: 



n-1 

an = a- aqnTn-i + a (^i^i-i - Qj^j-i) 

j=-oo 



-oo 



= a- aqnTn-i - a ^ {qj^j-i - Qj^j-i), 

j=n+l 
n-1 

j3n = (3- brnQn-i + ^ (^j^i-1 - rjQj-i) 

j = -QO 

+ 00 

= /? - &f„g„_i - ^ (^i^j-i - 

j=n+l 



n— 1 ^ — ~ +00 

C 



5 



1 — q',r,- -I- 1 — g^r,- 

j=— 00 J=n 



n-1 



d n i^^ = rfrr- 

J-l 1 - (7,r, J-1 1 



The equivalence of the two expressions for each auxiliary variable is re- 
lated to the fact that the three quantities X]j=^oo 'Y^i=-oo°''3li~^i 
rijL^oo ~ conserved as in the continuous-time case. These con- 

served quantities are assumed to be finite, and the last one should be nonzero. 
The substitution of these expressions into the first two equations in f l3.23p 
provides a global-in-space time discretization of the nonreduced Ablowitz- 
Ladik lattice (13.191) . In the continuum limit of time /i — 0, we obtain 

Q^n ^ a (lQnTn_ij (3n ^ (3 ^^nQ'n— 1; (^n ^ C, dji > d. 

Thus, with an additional but unessential constraint on the parameters, the 
time discretization fl3.23p reduces to the nonreduced Ablowitz-Ladik lattice 
(IXTUj) in the limit h ^ 0. 

The determinant of the 2x2 Lax matrix L„ fl3.20ap can be immediately 
computed as det L„ = 1 — qnTn- We compute the determinant of the 2x2 
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Lax matrix Vn (13.221) as 

det Vn{X) = X^h'^aCn + —h'^bdn 

A 

+ {ha - h'^a^qnTn-i + h^a(3n + /i^ag„_if„(i„ + hcn [1 + hon + ha (g„r„_i + g„r„_i)]} 
+ {hb- h%'^qn-irn + h%an + /i^6g„r„_ic„ + hdn [1 + h(3n + hh (g„-.ir„ + g„_ir„)] } 

+ h^ah + /i^c„(i„ + (1 + hqnTn^iCn + han) (1 + hgn-lVridn + h(3n) 

- h"^ {aqn - qnCn) {hvn - rndn) - h"^ {hqn-i - qn-idn) {arn-1 - r„_ic„) . 

Thus, the equahty fl2.10p combined with the boundary conditions (13.241) (or 
the streamhned version as stated above) imphes the set of relations 

ha - h^a^qnTn-i + h^ajdn + h'^aqn-ir^dn 

+ hcn [1 + han + ha (g„r„_i + g„r„_i)] = [ha + h'^a/3 + hc{l + ha)] A„, 
hb - h%\,_irn + h%an + h^qnT^-iCn 

+ hdn [1 + h(3n + hh (g„-ir„ + qn-iVn)] = [hb + h%a + hd{l + /i/?)] A„, 

h'^ah + h^Cndn + (1 + hqnVn-lCn + /ittn) (1 + hqn^-^lVndn + h(3n) 
- h'^ {aqn - qnCn) {hVn - rndn) - h"^ {hqn~l - qn-ldn) {oTn-l - Tn-lCn) 

= [h^ah + h^cd + (1 + ha){l + hl3)] A„. 

(3.25) 

Here, A„ is the quantity that satisfies the recurrence formula 

(1 - g„r„) A„ = (1 - g„r„)A„+i, (3.26) 

and has the normalized boundary value lim„^-|-ooA„ = 1. More precisely, 
only one of the two conditions lim„^_oo A„ = 1 and lim„_»_|_oo A„ = 1 is re- 
quired, in accordance with the boundary conditions for (q;„, /5„, c„, dn). Then, 
the other condition can be confirmed. Thus, A„ can be written explicitly as 




(3.27) 



Note that the first two equalities in (13.251) are valid even if a = or 6 = 
(cf. (13.231) ). The algebraic system (I3.25P essentially comprises three nontriv- 
ial relations for the three unknowns: «„, and A„. Before discussing the 
general case, we consider two special cases: a = 6 = 0orc = (i = 0. 
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• The case a = b = 0. 

In this case, the two auxihary variables, a„ and become n-independent, 
i.e., a„ = a, Pn = P- Thus, the discrete-time system (13.231) reduces to 

f I ^ 

^ (gn - gn) + (1 - g„r„)A„+i (cg„+i - dqn-i) - aqn + (3qn = 0, 

^ (jn - r„) + (1 - g„rn)A„+i {drn+i - cr„_i) - /?r„ + ar„ = 0, 

(1 - g„f„) A„ = (1 - g„r„)A„+i, hm A„ = 1 or hm A„ = 1, 

(3.28) 

and the algebraic system (13.251) simplifies to 

h^cd (1 - g„f„) (1 - g„_ir„_i) A^ - [(1 + ha) (1 + h(3) + h^cd 

- (1 + ha) hdqn-iVn - (1 + /i/3) /icg„r„_i] A„ + (1 + ha) (1 + /;,/3) = 0. 

(3.29) 

System (I3.28P is invariant under the following transformation: g„ i-^ f^^qn, 

^ fJ'~"'rn, he ^ fi^^hc, hd i— > fihd, ha i— > ha, h(3 t— > h(3, where is a nonzero 
constant. Using a similar transformation with respect to the time direction, 
it is possible to remove the parameters a and but we prefer to retain them. 
In terms of the normalized parameters 

he hd , , 

c:=q d:=- — , 3.30 

1 + ha l + h(3 ^ ^ 



(I3.29P can be rewritten as 

Cd (1 - qnTn) (1 - qn-lTn-l) A^ - 



1 + cd- dqn-iVn - cg„r„_i 



A„ + 1 = 0. 
(3.31) 



When cd ^ Q (and thus 7^ 0), this is a quadratic equation in A„, which 
has two solutions. The simplest way to reject the improper solution is to 
recall the asymptotic behavior of the Lax matrix Vn for small h (12. 8p . but 
here we take a different route. Let us first consider the "trivial" case where 
qn and r„ are zero for all n; if this is satisfied at some instant m = mo, then 
it holds true identically for any time m. Thus, the two solutions of (I3.3ip 
are given by A„ = 1, l/(c(i). The recurrence formula for A„ in (13.280 implies 
that these two solutions are unconnected, i.e., A„ takes the same value for 
all n. The solution A„ = l/{cd) can be discarded if 7^ 1, because it is 
inconsistent with the boundary condition for A„. Next, in the general case 
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where g„ and r„ are not identically zero, we assume that their amplitudes 
are always so small that their effects can be regarded as weak perturbations 
of the identically zero case. In particular, the value of A„ is restricted to a 
neighborhood of unity. Thus, we obtain the proper solution of the quadratic 
equation f l3.3ip as 

An = , ^ =, (3.32) 

1 + + V(1 + C„)'-4D„ 

with 

Cn '■= cd — dqn-iVn — C(7„r„_i, 
Dn := cd (1 - g„r„) (1 - g„_ir„_i) . 

Note that this solution is also valid for the linear case cd = 0. The decaying 
boundary conditions for g„ and r„ imply that lim„_>_oo A„ = lim„^+oo A„ = 1. 
When cd G M, the local expression fl3.32l) is valid only ii cd < 1. If > 1, 
the other solution of f l3.3ip should be adopted; alternatively, one can un- 
derstand the right-hand side of f l3.32p as being defined by a Taylor series 
for small c and d, and its analytic continuation. At the "threshold" value 
cd = 1, the discrete-time system (13.281) has the trivial solution g„ = dqn-i, 
= cr„_i, A„ = 1/(1 — qn-iTn-i). Thus, the discriminant of the quadratic 
equation fl3.3ip vanishes at cc? = 1, and the two solutions indeed intersect. 
Unless cci = 1, a unified expression for A„, 



A 



n 



l + Cn+{l-cd)M\ 



1 + Cn \ _ 4D„ 
i-cd ) (l-cd)2 



can resolve the sign problem of the square root. In addition, the above 
square root allows the Taylor expansion with respect to {qn^rn^qn-i^fn-i}- 
However, we do not use this unwieldy formula. 

The first and second equations in fl3.28p are linear in g„ and rVi, respec- 
tively. Thus, the forward time evolution that is explicit under the decaying 
boundary conditions can be expressed as follows (cf. f l3.30p ): 



l + ha f 1 + ha - 

-^—f^Qn + (1 - qnTn) K+l ( ~ ^ ^ /^/j ^^"+^ ^ ^""^ 

l + h(3 . . ^ ( l + hl3 

■r„ + (1 - g„r„) A„+i -— — — dr„+i + cr„_i 



+ ha \ 1 + ha 



(3.34) 



where A„ is given by fl3.32p with fl3.33p . In the simplest case of c = c = 
or d = d = 0, the forward time evolution f l3.34p is given by a simple rational 
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mapping. This fact was disclosed by Suris [7,19]. In the general case, using 
the new parameters c:=c[l + cd) and d:=d(^l + cd) , the mapping 
fl3.34p can be rewritten in a slightly simpler form: 



l + hp 
1 + ha 



ha 2(1- g„r 

qn + - 



Tn + 



Fn 

2 (1 - g„r„) 



\ + ha 
l + hl3 

l + h(3 
1 + ha 



1 ' 



rfr„+i + cr, 



n-l 5 



1 - dqnTn+l - CQn+lVn 



with 



+ Y (l - dqnTn+i - cqn+iTn) - 4:cd (1 - g„+ir„+i) (1 - g„r„) . 

That is, the first term cd of Cn+i is removed. When cd G M, this expression 
for Fn is valid only if — 1 < cd < 1; this corresponds to the range cd < 1/4. 
If cc? < —1 or > 1, the sign in front of the square root has to be changed. 

When h is real, (3 = a*, and d = c*, we can impose the complex con- 
jugacy reduction r„ = erg* with a real constant cr. In particular, setting 
a = —jS = —i'j/A'^, c = —d = —i/A"^, and r„ = —A'^q*, we obtain the fully 
discretized NLS equation 



h 



{q-n - Qn) + 



2 (1 + Z\2 



+1 



1 + C„+i' 



4D 



Z\2 



71+1 



{qn + Qn) = 0, 

(3.35) 



where 

Cn+l 
Dn+1 



(Z\2 



h^ lA^h _ lA^h ^ , 



/l2 



[l + A' \qn^,\') {l + A'\qn\'). 



The choice of 7 = 1 for the real parameter 7 is the most natural when taking 
the continuum limit, while the choice of 7 = simplifies the equation consid- 
erably. The aforementioned condition cd < 1 implies that (1 — 7^)/;,^ < 
this is automatically satisfied if 7^ > 1. At first glance, it is far from evident 
that the quantity in the square root is nonnegative: (1 + Cn+if — 4D„+i > 0. 
However, this inequality should hold true because the left-hand side repre- 
sents the discriminant of the quadratic equation in A„+i (cf. (13.311) ). In fact. 
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the reduction r„ = — z^^g* with > guarantees the auxihary variable A„ 
to be positive (cf. fl3.27p ). Thus, we also have the inequality 1 + Cn+i > 0. 
To summarize, qn+i and qn are not fully independent and satisfy the inequal- 
ity 1 + Cn+i > 1\JDn^\. Note that (13.351) with n ^ n + 1 determines 
from g„+2, and g„. 

Similarly to the continuous NLS equation, this full discretization is homo- 
geneous under the following weighting scheme: weight (Z\) = — 1, weight ih) = 
weight (g„) = 1. At present, it is unclear whether the corresponding one- 
parameter group of scaling symmetries, {A,h,qn) ^ {A/k,h/k'^,kqn), can 
define meaningful "self-similar" solutions to (13.351) . 

As mentioned previously, it might be aesthetically pleasing to eliminate 
the first "constant" term of C„+i by a suitable redefinition of the parameters. 
For example, in the simple case of 7 = and A'^ = 1, (I3.35P can be rewritten 
as 

i ,~ ^ , 2(1 + (g„+i + g„_i) 

T [qn - qn) + / , — = 0, (3.36) 

" 1 + i6Cn + \ (l + i6C„f - 46^Vn 



with 

Cn ■■= qnqn+l " Qn+iql, T^n ■= (l + \qn+lf) (l + \qnf) ■ 

Here, 6 := h/{l + h"^) is a new parameter. Because —l<h<l, the range 
of 5 is —1/2 < (5 < 1/2; at the end points 6 = ±1/2, the time evolution is 
trivial because cd = 1. Note that (13.361) is invariant under the transformation 
qn ^ (-l)"gn, S ^ -6. 

To express the backward time evolution explicitly, we only have to replace 
(1 — g„r„)A„_|_i in the first and second equalities of (13.281) with (1 — g„r„) A„, 
and then substitute the local expression (I3.32p with (I3.33p . 

In the complex conjugacy reduction wherein h is real, /? = «*, and d = c*, 
we can "normalize" the scaling of the dependent variable by setting r„ = — . 
In this case, we rewrite the first equation for g„ in (13.281) as 

(1 + ha*) qn - {1 + ha) + (1 + |gn|^)A„ {hcqn+i - hc*qn-i) = 0. (3.37) 

Moreover, (13.291) can be rewritten as 

|1 + haf + /i^|c|^ + (1 + ha) hc*qn-iq* + (1 + ha*) hcqnq*_i 

= |1 + haf + hW^ (1 + \qnf) (1 + A„. (3.38) 

Using (13.370 . we can replace 1/A„ and A„ in (I3.38P with rational expressions 
in qn as well as its shifts and complex conjugate. Thus, we obtain a rational 
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form of the fully discrete NLS equation, 

|1 + haf + h'^\c\^ + 2Re [(1 + ha) hc*qn-iq*] 
— h _i_ 7i |2 h -I- l~ P^l hcQn+i — hc*qn-i 

+ h^\c\' (1 + |g„_,n i^ + h^)ln-il + ha*)q. _ 

hcqn+i - hc*qn-i 

Surprisingly, this coincides with the double-discrete NLS equation proposed 
by Quispel, Nijhoff, Capel, and van der Linden [15, 17], up to a minor 
change of coordinates and parameters; despite its "elegance", this rational 
version has the drawback that the forward/backward time evolution cannot 
be uniquely determined. In addition, one cannot immediately recognize that 
fl3.39p reduces to the NLS equation in a continuous limit. Thus, we prefer 
our version, which is seemingly less elegant because of the existence of the 
square root but can define the unique time evolution properly and allow an 
easy-to-follow continuous limit. 

Actually, using a coordinate transformation, it is possible to "identify" 
our time discretization [ fl3.34p with d = c*, a = j3 = 0, and r„ = — g*] with 
the auto-Backlund transformation of the Ablowitz-Ladik hierarchy derived 
by Nijhoff, Quispel, and Capel [14]. However, their expression (see (19) 
in [14]) involves the indefinite sign ± in front of the square root and it is not 
clear how to understand and determine it. 

Once the local expressions for the auxiliary variables have been derived, 
we can normalize the 2x2 matrix Lax pair, Ln and Vn, so that det L„ and 
det Vn become equal to 1. Indeed, this is easily achieved by dividing Ln and 
Vn by V det Ln and ^/detVn, respectively (cf. (12.61) and (12.101) ). However, 
the normalized Lax pair not involving the auxiliary variables appears to be 
rather cumbersome and we do not present it here. 

• The case c = d = 0. 

In this case, the two auxiliary variables c„ and dn vanish. Thus, the 
discrete-time system (I3.23P reduces to 

^ {qn - qn) - aqn+i + hqn-i - a{qnrn-i - qn+irn)qn - anqn + qnf^n = 0, 

^ (jn - r„) - bVn+l + aVn-l - WnQn-l - Tn+iqnYn " (^nTn + ^nCtn = 0, 

an+i - an = a(g„r„_i - g„+ir„), 

(3.40) 
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where li m„^_ oo /5„) = (a, (3) or lim„^+oo ("n, Pn) = (a, The algebraic 
system (I3.25I) simphfies to 



1 — haqnTn-i + h(3n = (1 + h(3) A„ if a 7^ 0, 
1 - hbqn-irn + hun = (1 + ha) A„ if 6 7^ 0, 

hj^ah (1 - QnTn - qn-lTn-l) + (1 + hUn) (1 + h(3n 

= [h'^ab + (1 + ha){l + hp)] A„. 



(3.41) 



in; 



System f l3.40p is invariant under the following transformation: g„ i-^ /i"?n 
fn ^ f^~"fn, ha I— > fi'^ha, hb i-^ /la ^-^ ha, hf3 ^ hf3, where /z is a nonzero 
constant. 

It can be easily verified that the first equality in f l3.4ip is also valid for 
a = 0, and the second equality is valid for 6 = 0. Using the first and second 
equahties, the third equality in (13.411) results in a quadratic equation in A„, 

(1 + ha) (1 + h(3) Al - [{1 + ha) (1 + h(3) + h^ab - (1 + ha) haq-aTn-i 
- (1 + h(3) hbqn^iTn] An + h^ab (1 - qnVn) (1 - qn^l^n^l) = 0. (3.42) 

In terms of the normalized parameters 

ha f hb , ^ 

«:=^ b-=-, ^, 3.43 

l + hp 1 + ha ^ ' 

(I3.42P can be rewritten as 

Al - (^1 + ab- ag„r„,„i - 6g„_ir„ j A„ + ab{l- q^Vn) (1 - g„-ir„_i) = 0. 

(3.44) 

We recall that the Lax matrix Vn given by (13.221) with c„ = c?„ = is required 
to allow the asymptotic expansion (12. 8p for small h. This implies that the 
auxiliary variables an and determined by (I3.4ip are at most of order 0(1) 
and do not involve 1/h. Thus, one of the two solutions of the quadratic 
equation (I3.44p is rejected, and we obtain its proper solution as 



1 + A„ + J(l + Anf - 

An = , (3.45) 

with 

An := ab - ag„r„_i - bqn-iVn, 

_ _ (3.4o) 
Bn ■■= ab (1 - g„r„) (1 - g„-^ir„_i) . 
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The decaying boundary conditions for g„ and r„ imply that hm„^_oo A„ = hm„_^+oo A„ 
and consequently, lim„^_oo {an, Pn) = ^^^n^+oo (^n, Pn) = P)- When db G M, 
the local expression f l3.45p is valid only if afe < 1. If afe > 1, the other solu- 
tion of f l3.44p should be employed; alternatively, one can understand the 
right-hand side of f l3.45p as being defined by a Taylor series for small a and b, 
and its analytic continuation. At the "threshold" value db = 1, the discrete- 
time system (13.401) has the trivial solution g„ = ag„+i, r„ = 6r„+i, a„ = a, 
Pn = P- Thus, the discriminant of the quadratic equation fl3.44p vanishes 
at db = 1, and the two solutions indeed intersect. Unless db = 1, a unified 
expression for A„, 



l+^„H-(l-a£)^(^) 



' 5 



can resolve the sign problem of the square root. In addition, this expression 
allows the Taylor expansion with respect to {g^, ^n, Qn-i, ^n-i}- However, we 
avoid the use of this unwieldy form. 

Using the first two relations in fl3.4ip . we can rewrite a„ and Pn in the 
first two relations in f l3.40p in terms of A„, 



^ {in - Qn) + PQn " OiQu + + Ha-l) = 0, 



'■n 



^ {rn - Tn) + arn - Pvn + {-bvn+i + ar„_i) = 0. 



(3.47) 



This system is linear in g„ and r„. Thus, the explicit forward time evolution 
under the decaying boundary conditions can be expressed as follows (cf. 
dHHD): 



1 + ha 1 — g„r„ / 1 + /la - 

■.qn H 7 aqn+i - , ,^ ogr»-i 



l + hp l-g„r„ A l + 

'^n = ^ . rn H 1 I br n+i - , , , avn-i I , 



1 + ha An V 1 + /la 



(3.48) 



where A„ is given by (13.450 with (13.460 . 

In the simplest case ofa = a = 0or6 = 6 = 0, the quadratic equation 
(I3.44P in A„ reduces to a linear equation, and the forward time evolution 
(I3.48P is given by a simple rational mapping. This fact was discovered by 
Suris [7, 19]; actually, slightly prior to his work, relevant results in a prelimi- 
nary form were reported by Pempinelli, Boiti, and Leon [40]. 

When h is real, /? = «*, and b = a*, we can impose the complex con- 
jugacy reduction r„ = erg* with a real constant a. In particular, setting 
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a = —[3 = —i'y/A'^, a = —b = i/A^, and r„ = — Z^^g*, we obtain the fully dis- 
cretized NLS equation 

T- (gn - gn) + / , = 7^ + 5-) = O' 

(3.49) 

where 

iA^h _ iA^h ^ ^ 

Thus, this is equivalent to the fully discrete NLS equation (13.351) in the 
case a = 6 = 0, up to the space reflection n — >• — n; this correspondence can 
be readily noticed by considering the dispersion relation [3]. The choice of 
7 = 1 for the real parameter 7 is the most natural for taking the continuum 
limit, while the choice of 7 = simplifies the equation considerably. The 
aforementioned condition ab < 1 implies that (1 — 7^)/i^ < A^; this is auto- 
matically satisfied if 7^ > 1. Note that Qn-i and g„ are correlated and always 
satisfy the inequality 1 + v4„ > 2a/B^ if > 0. 

To express the backward time evolution explicitly, we have to rewrite 
A„ in (13.471) in terms of A„+i using (I3.26p . and then substitute the local 
expression (13.451) with (13.461) . By this procedure, the fully discrete NLS 
equation (I3.49P in the case c = d = can also be identified with the fully 
discrete NLS equation (13.351) in the case a = b = through the time reflection 
m —m and the sign inversion of h, h —>■ —h [3, 7, 19]. 

Note that the auxiliary variables appearing in the stage of time discretiza- 
tion can be interpreted, after a minor redefinition, as giving the fluxes of the 
first few conservation laws. Indeed, (I3.26P and (I3.40p provide the following 
conservation laws: 

A+ log (1 - g„r„) = A+ logA„, (3.50a) 
{QnTn-i) = A+ [(«„ -«)/« + g„r-„-i] , (3.50b) 
A+ (g„_ir„) = A+ [(/5„ -(3)/b + g„_ir„] . (3.50c) 

The cases a = and 6 = can be understood as limiting cases; the final local 
expressions for the fluxes do not contain 1/a and 1/6, and also hold true for 
these cases. Using ( Km . (Km . ([333]), ([335]), and ([33H]), we can express 
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the fluxes in fl3.50b|) and fl3.50cl) explicitly as 

1 + ha 



a 



ha 



A„ - 1 + hqn-lTn 



+ qnTn-l 



1 + ha 
ha 

1 + ha 
l + h^ 

h 

l + hp 



l + An-\(l + An)-ABn 



+ ab — aqnT, 



n-l 



QnXn—l 



2b (1 - g„r„) (1 - g„_ir„„i) 



+ b- g„r, 



n~l 



l + A„+^(l + A„)^-45„ 
2 (1 - g„_if„_i) (-6 + ag„f„_2) 



+ qnr n-l 



+ b 



Y {Pn - P) + qn-lTn = ^ [A„ - 1 + aqnTn-l] + g„-ir„ 

b hb 



l + h/3 
hb 

l + h(3 
1 + ha 

h 

1 + ha 



l+A^-Jil+AnY-ABr. 



2a (1 - g„r„) (1 - g„_ir„_i) 



ab - bqn-iTr 



+ qn-lTr. 



+ a - qn-lTn 



l + An + ^J{l + AnY-^Bn 
2 (1 - g„_if„_i) (-a + bqn^2rn) 



+ qn-lTr. 



+ a 



where A^ and B^ are given by (13.461) . Note that using (13.261) and (I3.48p . we 
can further rewrite these fluxes in many different forms. 

• The general case (a, b) ^ Q and (c, d) ^ 0. 

Let us consider the (nonreduced) discrete-time Ablowitz-Ladik lattice 
(I3.23P with (13.241) in the general case of (a, b) ^ and (c, d) ^ 0. Before 
dealing with the general case of the parameters a, 6, c, d, a, and /5, we first 
consider the special case where the conditions c = —a{^ 0), = — 6(7^ 0), 
and (3 = —a are satisfied. This is a particularly interesting case from the 
perspective of applications such as numerical experiments, because the linear 
part of the equations for g„ and r„ has symmetry with respect to the time 
reflection m —m combined with the sign inversion of h. Moreover, the 
symmetry with respect to the space reflection n —n can also be achieved 
by imposing the further condition b = —a [1,3,5,19,41]. 
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In the case of c = —0(7^ 0), d = — 6(7^ 0), and /5 = —a, the algebraic sys- 
tem fl3.25l) can be reformulated as 

(1 - hhqn-irnAn + h[5n) - A„ (1 - /iag„r„_iA„ + han) 

= haqnTn-i + [-2ha + ha {qnTn-i + qnTn-i)] A„ + haq^r^-iAl, 

(1 - haqnTn-iAn + /ia„) - An (1 - hbqn-iVnAn + hPn) 

< = hbqn-iVn + [2ha + hb (g„_ir„ + g„-ir„)] A„ + /i6g„_if„A^, 

(1 - haq^Tn-iAn + /ia„) (1 - hbq^-irnAn + h(3n) 
= -h^ab (1 - g„r„ - g„_ir„„i) - h'^ab (1 - g„r„ - g„_ir„_i) A^ 
+ [/i^a6 (2 + g„r„ + g„r„ + §n_ir„_i + g„_if„_i) + (1 + - ha)] A„, 

(3.51) 

with Cn = —aAn and (i„ = —bAn- Using the first two relations, we can express 
1 — haqnfn-iAn + han and 1 — hbqn-iVnAn + hj3n in terms of A„. Substitut- 
ing these expressions into the third relation, we obtain a sixth-degree equation 
in only A„. However, this sextic equation turns out to be too complicated 
for exact treatment. Thus, we need to take a different route. One sensible 
approach is to solve the algebraic system (13.511) approximately; for exam- 
ple, using perturbative expansions with respect to /i, we obtain approximate 
expressions for the auxiliary variables as follows: 

< ft 

an - ag„r„_iA„ = « - - (g„ + g„) (r„_i + r„,_i) + 0{h), 

b ^ 

< Pn - bqn^iTnAn = -a - - (g„_i + g„„i) (r„ + r„) + 0{h), 

An = l-^{qn + qn) (^„-l + Tn-l) - ^ {qn-1 + ^n-l) (^n + r„) + 0(/l^). 

(3.52) 

Combining (I3.52p with (13.231) for c = —a, d = —b, and P = —a, we obtain 
a high-precision numerical scheme that is "almost integrable" (cf. [3,5,41]). 
We can also consider the next order in h and include higher-order corrections. 

We now discuss the general case where no particular conditions on the 
parameters a, b, c, d, a, and (3 are assumed. Instead of exploiting the compli- 
cated algebraic system (I3.25p . we explain how to recreate the discrete-time 
Ablowitz-Ladik lattice (I3.23P by the composition of a few more elementary 
time evolutions such as those studied earlier for the special cases a = 6 = or 
c = d = 0. To this end, we recall some basic properties of integrable hierar- 
chies. That is, under suitable boundary conditions, an integrable hierarchy 
generally comprises infinitely many commuting flows. Any composition of 
these flows also belongs to the same hierarchy; each of the flows is uniquely 
determined by specifying its linear part, or equivalently, its dispersion rela- 
tion. Under the boundary conditions (I3.24p or, equivalently, the streamlined 
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version, the discrete-time Ablowitz-Ladik lattice fl3.23p has the following lin- 
ear leading part: 



J- {qn - qn) + cqn+1 + /5g„ + hqn-i - ag„+i - aqn - ~ 0, 

^ (3.53) 

T {rn - Tn) + dvu+i + ar„ + OTn-i - hr^+i - (3rn - cr„_i ~ 0. 
h 

Substituting the ansatz g.„ ~ A^", g„ ~ X'^^uj, r„ ~ A"^", r„ ~ A"^"a;~^ into 
fl3.53p . we obtain the corresponding dispersion relation [3] 



Note that the Lax matrix Vn f l3.22l) approaches a diagonal matrix as n ^ ±oo 
(cf. (13 ■241) ). and oj (A^) is equal to the limit value of the ratio between the 
diagonal elements [3] ; this is natural because the n-independent overall factor 
of Vn plays no role in the zero- curvature condition (12. 6p . 

For notational convenience, we express the forward time evolution defined 
by (I3.23P with (I3.24p as the mapping 

/(a, h,c,d,a + \,I3 + \) : (g^, r„) (g„, f„) . 

Note that only the ratio of the six arguments in / matters in identifying the 
corresponding flow as well as its linear part. Let us discuss how to factorize 
this mapping into a composition of simpler mappings belonging to the same 
class. A rigorous and persuasive discussion can be given on the basis of the 
scattering data used in the inverse scattering formalism. The scattering data 
are space-independent quantities that have very simple time dependences de- 
termined by the dispersion relation, typically, -R(A) = oj (A^) -R(A) [3]. Thus, 
their time evolution is essentially linear and allows a commutative superpo- 
sition of different flows. Conversely, the decomposition of a complex flow 
into elementary flows can be uncovered by examining the factorization of the 
dispersion relation, that is, whether oj (A^) can be factorized into a product 
of simpler quantities such as Ui (A^) 0J2 (A^) or uji (A^) UJ2 (A^) u^, (A^). If the 
reader is unfamiliar with the inverse scattering method, the same guideline 
can be derived by noting the one-to-one correspondence between (the linear 
part of) a flow and its dispersion relation. One can also understand this 
guideline "honestly" at the level of the associated linear problem (cf. the 
second equation in (12.50 and the sentence beneath (13.540 ): that is, by consid- 
ering whether the Lax matrix Vn can be factorized into a product of simpler 
matrices such as vP (A) V^^^ (A) or V^^^ (A) vP (A) V^^^ (A) [1,6,7,19]. 
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By elementary computations, the right-hand side of fl3.54p can be factor- 
ized into the product of two or three simpler quantities (as functions of A) 
as follows: 

1 + h (aA^ + a + ^) 
1 + h (cA2 + f3 + ^) 

^ 2hd/X'^ 

1 + ha + + hay — AK^ad + 2ha\^ i+ha+^{i+ha)'2-ih'^ad 



X 



l + h(5+ v/(l + hPY - 4h%c + ^ 1 



2/icA2 



l+h/3+y/{l+hl3)'^-4:h^bc 

(3.55a) 



1 I 2haX'^ 

i+ha+y/{i+ha)2-4h^ad 1 + ha + a/(1 + ha^ - Ah^ad + ^ 

X 



^^^'^^ l + hj3+ J(l + hpy -Ah%c + 2/icA2 

l+hl3+y/{l+hl3)^-Ah'2bc ^ 



(3.55b) 



2 _| 2ha\'^ 

1 + ha + a/(1 + hay — AK^ad l+ha+yj{l+ha)'^~ih'^ad 

1 + + vTl + h(3Y - Ah^W ^ 1 + J^^M 



l+h/3+y/ (l+hfSy-'-ih'^bc 



1 + 



2hd/X^ 



^ 2 I 2hcx^ • (3-55c) 

l+h|3+^y{l+h|3)^-4h^2bc 

Thus, the aforementioned guideline implies the corresponding decomposition 
of the mapping /(a, b,c,d,a + j^, P + j^) into simpler ones, i.e., 

/(a, 6, c,d,a + j^,(3+l) 

= f{a, b, 0, 0, a' + i, /?' + i) o /(O, 0, c', d', i, i) (3.56a) 

= /(a', 6', 0, 0, i, i) o /(O, 0, c, d, a' + ^ + D (3.56b) 

= /(O, 0, 0, 0, a' + lP' + l)o f{a', b', 0, 0, i, i) o /(O, 0, c', rf', I i), 

(3.56c) 



where 



a' 



c' 



2a 


l + ha + a/(1 + hay 


-Ah^ad 


2c 




l + v/(l + /i/9)2 


-Ah%c 


v/(l + hay - Ah^ad - 


- 1 + ha 


2h 



d' ■-- 



1 + + ^{l + hpy -AK^bc 
2d 



l + ha+ a/(1 + hay - Ah^ad' 



. , ^, v/(l + hf3y - Ah%c -l + h(3 



2h 
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Note that the order of composition in fl3.56p does not matter, as the order 
of muhiphcation in (13.551) does not matter. Interestingly, the time evolu- 
tion defined on the six points {n + i,m + j), i = —1, 0, 1 and j = 0, 1, is now 
written as the composition of two time evolutions, each of which is a previ- 
ously described four-point scheme, and, if any (cf. f l3.56cl) ). one trivial time 
evolution. Moreover, if h is real, b = a*, d = c*, and f3 = a*, the original 
mapping /(a, b,c,d,a + j^, P + ^) allows the reduction of the complex conju- 
gate r„ = erg* (cr: real), and this reduction is consistent with every component 
mapping involved in (13.561) . 

Using the dispersion relation (I3.54p . we can readily find the inverse map- 
ping of f {a, b,c,d, a + j^, P + j^). Indeed, the trivial identity 

1 + h {aX^ + a + ^ 1 + h (cA^ + P + ^) _ ^ 
1 + h (cA2 + (3 + ^) 1 + h (aA2 + a + ^)~ 

implies the nontrivial formula 

/-^(a, 6, c, c?, a + i, /3 + i) = /(c, d, a,b, (3 + ^, a + I). 

In particular, this confirms the already mentioned correspondence between 
the case a = b = and the case c = d = through the time reflection. 

Actually, the right-hand side of (13.541) is already the product of the nu- 
merator and the inverse of the denominator; this implies the decomposition 

/(a, 6, c, a + i, /? + i) = /(a, 0, 0,d,a + \, o /(O, 6, c, 0, i, /5 + 

= /(O, 6, c, 0, /? + i) o /(a, 0, 0, c?, « + i, i). 

Thus, using this formula, one can further decompose (I3.56P into "elementary" 
flows (cf. [1,7,19]); however, such flows do not maintain the complex conju- 
gacy relation between g„ and r„, which may become a serious bottleneck in 
practical applications. 

3.3 The Volterra lattice 

The Volterra lattice 

Un,t = UniUn+l - Un-l) (3.57) 

is obtained from the Ablowitz-Ladik lattice (I3.19P through the reduction 
a = 6 = 1, g„ = u„ — 1, and r„ = —1. Thus, the time discretization of the 
Volterra lattice can be obtained from the discrete-time Ablowitz-Ladik lat- 
tice in subsection 13.21 through the reduction g„ = m„ — 1 and r„ = — 1 to- 
gether with a = b, an = Pn+i, and Cn = dn- The discrepancy between this 
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reduction and the decaying boundary conditions (13.241) is nonessential in this 
regard. However, in the following, we prefer to consider the Volterra lattice 
independently as an interesting example. The Lax pair for the continuous- 
time Volterra lattice (13.571) is given by [1,42] 



A 



Un 





n Alt„ 



-A Un-l 



(3.58a) 



(3.58b) 



Indeed, the substitution of (I3.58P into the zero-curvature condition (12. 2p 
results in (I3.57p . 

In the discrete-time case, we assume the Lax matrix Vn of the following 
form: 



Vn = I + h 



A^a + a„ + Unh„ 
-A (a + hn) 



A (aUn + Unbn) 
-X%n + Oin~l + Un-lbn-1 



(3.59) 



Here, an and 6„ are auxiliary variables. The zero-curvature condition (12. 6p 
for the Lax pair (I3.58ap and (I3.59P is equivalent to the following system of 
partial difference equations (cf. [5]): 



^ {Un - Un) = an+lUn - Unttn-l + Un+l&n+lMn " M„Mn-l&n-l, 

an+l -an = a {Un+l - Un) , 
Unbn = bn+lUn- 



(3.60) 



The general form (I3.60p of the time-discretized Volterra lattice is integrable 
for matrix-valued dependent variables, but in the following, we consider only 
the case of scalar dependent variables. In view of (12. 9p . we impose the fol- 
lowing boundary conditions for Un, ctn, and 6„: 



lim Ur. 

n— >itoo 



lim an 

n— >±oo 



0, 



lim b„ 

n— >itoo 



b. 



(3.61) 



The boundary value of m„ is normalized by scaling a and b, and the boundary 
value of On is set as zero by redefining h. To be precise, the boundary 
conditions (I3.6ip contain redundant information. Indeed, it can be shown 
that the auxiliary variables have the same limit values for n — oo and 



n 



+00. Therefore, it is sufficient to assume either lim^ 



(0,6) 



or lim„ 



aUr 



reduces to 



(a„, 6„) = (0, b). In the continuum limit of time /i ^ 0, we obtain 
a and 6„ b. Thus, in this limit, the time discretization (I3.60p 

Un,t = (a + b)Un{Un+l - Un-l), 
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which is, if a + 6 7^ 0, equivalent to the Volterra lattice (13.571) . up to a redef- 
inition of t. Note that in the case a + b = 0, (I3.60p has the trivial solution 

Un = Un, an = aUn - a, bn = b. 

The determinant of the 2x2 Lax matrix L„ f l3.58ap can be immediately 
computed as detL„ = Un- Using the recurrence formulas for a„ and 6„ in 
(I3.60p . we can rewrite the Lax matrix Vn (13.590 as an ultralocal form with 
respect to the auxiliary variables a„ and bn- Thus, its 2 x 2 determinant is 
computed as 

detK(A) 

= [l + h (A^a + a„ + [l + h (-A^6„ + a„ - aUn + aUn-i + &„m„-i)] 

+ /i^A^ {aUn + Unbn) (a + bn) 
= -X'^h'^abn + X^h [{a - bn) (1 + /ia„) + ha?Un~i + habn + Un + Un)] 

+ [1 + /l (a„ + Unbn)] [l + h [an - aUn + aUn-l + bnUn-l)] ■ 

Therefore, the equality fl2.10p combined with the boundary conditions fl3.6ip 
(or the streamlined version as stated above) implies the set of relations 

' bn = bAn, 

{a - bn) (1 + hOn) + ha^Un-l + habn {Un-l +Un + Un) 

= (a — 6 + ha^ + 3/ia6) A„, 

[1 + /l (a„ + Unbn)] [l + h {an — CLUn + CLUn-1 + bnUn-l)] = {l + hb)'An. 

(3.62) 

Here, is the quantity that satisfies the recurrence formula 

M„A„ = M„A„4.i, (3.63) 

and has the normalized boundary value limn^j^oo-^n = 1- More precisely, 
only one of the two conditions lim„^_oo A„ = 1 and lim^^+oo = 1 is re- 
quired, in accordance with the boundary conditions for a„ and 6„. Then, the 
other condition can be confirmed. Thus, A„ can be written explicitly (but 
globally) as 

n—l ~ +00 

An= n -=n- 

^■^ Uj u,- 

jr'=— oo j='n 

Note that the first equality in (13.620 is valid even if a = (cf. (13.600 ). The 
algebraic system (13.620 essentially comprises two nontrivial relations for the 
two unknowns and A„. 

Because of the nonzero boundary value of m„ (I3.6ip . it is rather cum- 
bersome to identify the linear leading part in (I3.6O0 and the corresponding 
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dispersion relation (cf. subsection I3.2p . Instead, we can compute the asymp- 
totic form of the Lax matrix Vn fl3.59p asn^— ooorn^+oo and consider 
its factorization (cf. [1,6,7,19]). The results imply that the cases a = and 
6 = form the basis for the general case; that is, the time evolution for gen- 
eral a and b is equivalent to the order-independent composition of two time 
evolutions corresponding to a = and 6 = 0, respectively. In addition, these 
two cases are related to each other through the space/time reflection, as we 
will see below. 

Now, we consider the two fundamental cases: a = or 6 = 0. 

• The case a = 0, 6 7^ 0. 

In this case, the auxiliary variable a„ vanishes. Thus, the discrete-time 
system (13.601) reduces to 

Un) = hkn+lUn (m„+1 - M„-l) , 

(3.64) 

u„A„+i, lim A„ = 1 or lim A„ = 1. 

71— » — 00 71— > + 00 

Note that in the case hh = —1, (I3.64p has the trivial solution Un = Un-i, 
An = 1/un-i- The algebraic system (13.620 simplifies to a quadratic equation 
in A„, 

(1 + hmnkn) (1 + hhUn-lK) = (1 + hhf An, 

or equivalently, 

{hhfUnUn-lAl - [(1 + hhf - hh {Un + M„-l)] A, + 1 = 0. (3.65) 

The asymptotic behavior (12.80 of the Lax matrix Vn implies that the proper 
solution of (I3.65P is given by 

(1 + hhf An = , ^ ^ ^ ^ --, 

1 - e (m„ + Un-l) + e{Un + Un-l)f - 4:t'^UnUn-l 

(3.66) 

where e := hh/ (1 + hhf. Thus, if hh G M, then e < 1/4. In this case, the 
local expression (13.660 is valid only if — 1 < /ifo < 1, which covers the range 
e < 1/4. The case hh = 1, corresponding to e = 1/4, involves a subtle sign 
problem of the square root, which we do not discuss here. If {hhf > 1, 
(I3.66P is inconsistent with the boundary conditions, and the other solution of 
(I3.65P should be adopted. In any case, the boundary conditions for m„ imply 
that lim„^_ooA„ = lim„^+ooA„ = 1. As hh approaches 1, more restrictive 
conditions have to be imposed on the u„ to ensure a consistent choice of the 
solution of (I3.65p . 
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Substituting (I3.66P into the first equation in fl3.64p . we obtain a time 
discretization of tlie Volterra lattice fl3.57p in tlie local form, 



~e " 1 / ~ 2 ~ ' ^ 



This is an explicit scheme under the boundary conditions lim„^-|-ooW„ = 1, 
because the value of is uniquely determined by «„, Un+i, and Un-i- Note 
that (13.670 can be rewritten as 

Un — Un 2 (Un+lUn — UnUn-l) /„ 

e ^ 1 / _ 2 ~ " 

1 + e - m„) + y [1 + e - - 4eM„+i 

If /i6 G M and the u„ are nonzero and real-valued at the initial time, then 
fl3.64p implies that the auxiliary variable A„ is always real-valued. Con- 
sequently, the discriminant of the quadratic equation (13.651) must be non- 
negative. Thus, if e < 1/4, M„ approaches 1 sufficiently smoothly and fast 
as n — s> ±oo, and the are close to 1 at the initial time, then the real- 
valuedness of m„ is preserved under the time evolution of the discrete-time 
Volterra lattice fIXUTI) or (KEE\\ . 

To express the backward time evolution explicitly, we only have to replace 
An+iUn in the first equation of (I3.64p with A„m„ and then substitute the local 



:. (3.69) 



expression fl3.66p . The resulting equation is 

Un - Un 2Un (^n+l " ^n-l) 



1- e{Un + Un-l) + y[l -e{Un + M„-l)] - Ae'^UnUn-l 

• The case 6 = 0, a 7^ 0. 

In this case, the auxiliary variable bn vanishes. Thus, the discrete-time 
system (13.601) reduces to 

T (Un - Un) + (^n - Un) = a (u„+iU„ - UnUn^l) , 

h ^ (3.70) 

— = a (wn+i — Un) , liui = or lim = 0. 

n— >— oo n— >+oo 

Note that in the case ha = —1, (13.701) has the trivial solution 1 + han = Un+i, 
Un = Un+2- Thus, we can assume ha ^ —1. The algebraic system (I3.62p 
simplifies to 

1 + han + haun-i = (1 + ha) An, 
(1 + han) (1 + han — haun + haUn-i) = A„. 
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By eliminating A„, (13. 711) reduces to a quadratic equation in 1 + 
(1 + ha) (1 + han)'^ — [1 + /la (1 + ha) {un — Un~i)] (1 + han) — haun^i = 0. 
The proper solution of this quadratic equation is given by 



1 + F„ - G„ + J{1 + F„- G.af + 4G„ 

1 + = j^tth ■ 

with 

Fn := /ia (1 + ha) Un, 
Gn := ha (1 + /la) Un~i- 

When ha G M, the local expression (13.721) is valid only if ha > —1/2; the 
borderline case ha = —1/2 is excluded from our consideration, li ha < —1/2, 
the other solution of the quadratic equation should be employed. Unless 
ha = —1/2, a unified expression for 1 + 



^ , _ l+Fn-gn + (l + 2Mv/(^-^) +(Tji^ 

+ " 2(1 + H 

can resolve the sign problem of the square root, but it is unwieldy and 
looks unattractive. Note that the boundary conditions for m„ imply that 
hm„^_oo an = lim„^+oo "n = 0, and consequently, lim„^_oo An = lim„^+oo An 

Substituting the local expression (13.721) for 1 + han into the first equa- 
tion of (13.701) . we obtain a time discretization of the Volterra lattice (13.571) . 
In terms of the new parameter 6 := ha {1 + ha), we can write this time dis- 
cretization as 

Un ~ Un 2 (Un+lUn — UnUn-l) 



l + 5{Un- Un-l) + + 5 {Un ~ Un~l)f + ^5Un~l 



(3.73) 



Under the boundary conditions lim„^±oo ^n = 1, this is an explicit scheme 
in the forward time direction. Note that if ha > —1/2, then 6 > —1/4. If 
this inequality is satisfied, Un approaches 1 sufficiently smoothly and fast 
as n ^ ±oo, and the m„ are close to 1 at the initial time, then the real- 
valuedness of u„ is preserved under the time evolution (I3.73p . 

To obtain the backward time evolution explicitly, we rewrite an in the 
first equation of (I3.70p as ttn+i — « (wn+i — Un) and then substitute the local 
expression (I3.72p . The resulting equation is 

Un-Un 2 [Un+lUn - UnUn-l) 



S (Un+l - Un) + 5 (Un+l - Un)f + 4:SUr. 
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(3.74) 



It is easy to see that fl3.73p is equivalent to fl3.68l) (or (13.691) ) through the 
space (or time) reflection and the identification 6 ^ — e. In the same man- 
ner, fl3.74p can be identified with f l3.68p through the time reflection. Thus, 
the forward/backward time evolution in the case 6 = corresponds to the 
backward/forward time evolution in the case a = 0, up to a redefinition of 
the parameters. 



Ultradiscretization [12,13] . We propose an ultradiscrete analogue of the time- 
discretized Volterra lattice in the case 6 = 0. For the forward time evolution, 
we first rewrite (13.730 as 



Ur. 



1 + 6{Un - Un-l) + 26Un+l + \/ 1 + 2(5 («„ + M„_l) + 5^ ( 



l + 6{Un + Un-l) + Y 1 + 25 (U„ + Un-i) + 5^ (m„ - U^-lf 

(3.75) 

Then, we assume 5 > so that the positivity of the dependent variable, 
Un > 0, Vn G Z, can be preserved in the time evolution. By reparametrizing 
the parameter and the dependent variable as 6 =: exp{—L/e) and m„ =: exp([/„/e), 
respectively, and taking the logarithm, the above equation becomes 



Un~Un=e log 



-elog 
Xn := l + 2e"^ +2e 



1 + e^^ - e s + 2e" 



1 + e^ + e^ + JXn 



+ e £ — e £ 



(3.76a) 
(3.76b) 



For brevity, the e-dependence of the [/„ is suppressed. When taking the 
limit e —>■ +0, the well-known formula [13,43,44], 



> e £ = max (Aj) , 



where Aj{e) are real functions allowing the one-sided limit lim^^+o Aj{e) =: Aj, 
does not apply directly; however, the core idea behind this formula is still 
valid. After simple consideration of the cases Un I Un-l when [/„ > L, flXTHj) 
is shown to reduce, in the limit e +0, to 



Un-Un = max f{Un, Un-l), Un+1 - L] - max [0,Un- L, Un- 



L]. 



(3.77) 
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Here, Un '■= linig^+o ?/„(£:) is a real variable, and the boundary conditions 
lim„^-i-oo w„ = 1 translate into lim„^-|-oo Un = 0. The function Un-i) is 

defined by 



fiUn, U, 



n—1) 



lim eloEf 



Un-L 

1 + e^^ 



Un - L, Un> Un-1 and Un> L 

Un-1 - Un-1 + g{Un-2, Un-l) , Un = Un-1 > L 
0, Un < Un-1 OI Un< L 



>o. 



where g{Un-2, Un-i) is defined in a similar manner as 



g{Un-2, Un-i) := lim elog 



1 + e £ - e e + VX„_i 



Un-2 - L, Un-2 > Un-1 and Un-2 > L 

fiUn-2, Un-3), Un-2 = Un-1 > 

0, Un-2 < Un-1 or Un-2 < L 



> 0. 



One might naively think that f{Un, Un-i) at f/„ = Un-i > L is given by (f/„ — L)l2] 
however, this is not true in general, because the equality f/„(+0) = t/„_i(+0) 
does not imply Un{s) = < £ -C 1. As indicated above, the correct 

value of f{Un, Un-i) can be computed recursively using the formulas 

/(f/„+l, Un) = Un-Un + g{Un-U Un) , g{Un, Un+l) = f{Un, Un-l). (3.78) 

We now briefiy explain how these formulas can be derived in the most general 
case, i.e., without any assumptions on the arguments such as f/„+i = Un- 
vote that (13.781) can be regarded as a conservation law. Thus, in view of the 
zero boundary conditions for [/„, the following global expressions are valid if 
L > 0: 



oo oo 

fiUn+1, Un) = - = 

i=o j=i 

oo 



+2j — Un+2j j , 



g{Un, Un+l) — (^Un-2j-l — Un-2j-l^ — (^Un+2j-l — Un+2'j-l 

3=0 j=l 
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To verify fl3.78l) . we compare fl3.75p with fl3.63l) . wherein the expression for 
A„ is given by the first equation in (I3.7ip with (13. 721) . This results in the 
nontrivial relation 

l + 6{Un- Un+l) + \Jl + 25 {Un+1 + M„) + 5^ ( Un+l ^ Un) 



= 1 + 6{Un- Un-l) + y 1 + 25 («„ + Un-l) + 5^ (Wn " ^n-l) • 

It is an easy exercise to rewrite this relation by rationalizing its numerator 



as 



1 + 5 {Un+l - Un) + \Jl + 25 {Un+l + M„) + 5'^ ( 

Un+l ~ Un) 

= — 1+6 {Un-l - Un) + \ 1 + 26 {Un + Un-l) + 6^ (m„ - U„_l)' 
Un L 

By substituting 6 = exp{—L/e) and Un = exp{Un/s), and taking the loga- 
rithm, the above two equalities reduce, in the limit e +0, to the second 
and first equalities in f l3.78p . respectively. 

The parameter L may depend on the discrete time coordinate m G Z. The 
commutativity of the ultradiscrete flows defined by f l3.77p for distinct values 
of L remains to be established. In this connection, we conjecture that fl3.77p 
and its time reversal (see below) comprise an ultradiscrete analogue of the 
KdV hierarchy, wherein the parameter L labels each flow in the hierarchy. It 
is also beyond the scope of this paper to investigate the relationship between 
the ultradiscretized Volterra lattice f l3.77p and the "box and ball system" of 
Takahashi and Satsuma [45]. 

For the backward time evolution, we rewrite (13.740 as 



Un 
Un 



l + 6{Un- Un+l) + 26Un-l + J I + 26 {Un + Un+l) + 6^ ( Un Un+l) 



l + 6{Un + Un+l) + \jl + 26{Un + Un+l) + 6"^ (m„ - 



Thus, in the same way as that for the forward time evolution, we obtain an 
ultradiscrete analogue of (I3.74p . 

Un-Un = max (g{Un, Un+l), Un-l ~ ^) ~ - L, Un+l " ^) , 

which is related to (13.770 through the combined space and time reflection (cf. 
flX7i|) and (Km ). 
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3.4 The modified Volterra lattice 



The modified Volterra lattice 



'1 + ?n)(?n+l -Qn~l] 



(3.79) 



where g„ is a scalar dependent variable, was introduced by Hirota [46] (also 
see [4,47,48]); it is obtained from the Ablowitz-Ladik system (13.191) through 
the reduction a = b = 1 and r„ = —Qn- The overall coefficient of the nonlinear 
terms, including its sign, is nonessential at the level of the equation and the 
associated Lax pair, because it can be changed by rescaling g„. Note that 
(I3.79P is invariant under the transformation g„ i-^ (— l)"'g„, t i— >■ —t. The Lax 
pair for the continuous-time modified Volterra lattice (13.791) is given by 



A 



1 

X J 



QnQn-i Ag„ + fqn-1 



1 



(3.80a) 
(3.80b) 



Indeed, the substitution of fl3.80p into the zero-curvature condition (12. 2p 
results in (I3.79p . 

In the discrete-time case, we consider the reduction a = 6, a„ = Pn, Cn = dn 
and Tn = —q-n of the Lax matrix Vn (I3.22p . Thus, we obtain 



Vn = I + h 



X^a — qnQn-lCn 


X{aqn - QnCn) 


A(-ag„_i + CnQn-l) 
+ i(-ag„ + QnCn) 


A Cjj QnQn—l^^n 



(3.81) 



where a„ and c„ are the auxiliary variables. The zero- curvature condition 
(12. 6p for the Lax pair (I3.80ap and (I3.8ip amounts to the following system of 
partial difference equations: 



h 



{qn - Qn) - aqn+1 + aQn-l - On+l^n + ^nCtn 

+ g„+ic„+i(l + ql) - (1 + g„^)c„g„_i 



a {-qnqn-i + qn+iqn 



i + 9„)c„ = (i + g„)c 



n+1- 



(3.82) 



The discussion in subsection 13.21 implies that the general case of a 7^ and 
c„ 7^ is equivalent to the composition of two simpler cases: the case a = 
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and a„ = and the case c„ = 0. Moreover, the latter case can be related 
with the former case through the time reflection. Thus, in this subsection, 
we only consider the elementary case a = 0, a„ = 0, and c„ = — under the 
boundary conditions 

lim qn = 0, lim A„ = 1. (3.83) 

n— >itoo n— >itoo 

Note that the boundary value of c„ is set as —1 without any loss of generality; 
it can be changed to any nonzero value by rescaling the "step size" parameter 
h. To be precise, the boundary conditions (13.831) contain redundant informa- 
tion. Indeed, it can be shown that A„ has the same limit value for n —oo 
and n —>■ +oo. Therefore, it is sufficient to assume either lim„^_oo A„ = 1 
or lim„_^+ooA„ = 1. Now, the general system (I3.82p reduces to the simpler 
system 

T &n - qn) = (1 + ql)^n+l (Wl " ^n-l) , ,^ 

h (3.84) 

{l + q^)An= {l + ql)An+i. 

The second relation in (I3.84p implies that the global expressions for the 
auxiliary variable A„ in terms of the Qn are 



A. 



1 + Q? 11 1^ 



+ g/ ^J-l + q 



The substitution of each expression into the first relation in (I3.84p provides a 
global-in-space time discretization of the modified Volterra lattice (13.790 . If 
/i G M and the g„ are real- valued at the initial time, then the real-valuedness 
of Qn is preserved under the discrete-time evolution. 

Let us resolve this nonlocality. The determinant of the 2x2 Lax matrix 
Ln (I3.80ap is given by detL„ = 1 + g^, while the determinant of the Lax 
matrix Vn (I3.8ip with a = 0, a„ = 0, and Cn = — A„ is computed as 

detK(A) =- (A^ + l/A^) hAn 

+ h^l + q^) (1 + Al + 2/ig„g„_iA„ + L 

Thus, the equality (I2.10p combined with the boundary conditions (I3.83P (or 
the streamlined version as stated above) leads to the quadratic equation in 

An, 

(1 + q^) (1 + Al-{l + h^- 2/ig„g„_i) A„ + 1 = 0. (3.85) 
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By recalling the prescribed asymptotic behavior fl2.8p of the Lax matrix Vn 
for small h, we obtain the proper solution of this quadratic equation in A„, 



1 + /i2 - 2/ig„g„„i + J(l + h^- 2hqnqn-iY - ^h^ (1 + q^) (l + qLi) 



l + K^ - 2/ig„g„_i + y (1 - K^f - Ah (g„ + hq^-i) (hq^ + g„_i) 

(3.86) 

When /i G M, the local expression fl3.86p is valid only ii —1 < h < 1. li > 1, 
the other solution of (13.851) should be used. Unless h = ±1, a. unified expres- 
sion for A„, 

2 

A„. = 



l + h"^ - 2/ig„g„_i + {1 - h"^) ^ 1 - j^^^ {Qn + hq^-i) (/ig„ + g„_i) 

can resolve the sign problem of the square root, but we do not use this 
form. In any case, the decaying boundary conditions for g„ imply that 
lim„^_oo A„ = lim„^+oo A„ = 1. Thus, the boundary conditions for A„ given 
in fl3.83p are compatible. 

Substituting (I3.86P into the first equation in (I3.84p . we obtain a time 
discretization of the modified Volterra lattice (13.791) in the local form, 

^ (?n - qn) 

2(1+ ql) {qn+1 - qn-l) 



l + K^ - 2hqn+iqn + \/(l - h^f - Ah (g„+i + /ig„) (/ig„+i + g„ 



(3.87) 



The right-hand side of (I3.87P does not involve g„, and thus the forward time 
evolution is explicit under the decaying boundary conditions. Similarly to the 
continuous-time case, (I3.87P (or (I3.84p ) is invariant under the transformation 
qn ^ (— l)"''?n, h I— i> —h. \i —1 < h < 1, the real-valuedness of qn is preserved 
under the discrete-time evolution (I3.87P : the discriminant in the square root 
is nonnegative as long as we start with sufficiently small real-valued g„ at 
the initial time. To express the backward time evolution explicitly, we only 
have to replace (1 -f- g^)A„+i in the first equation of (I3.84p with (1 + q^) A„ 
(cf. the second equation) and then substitute the local expression (I3.86p . 

In terms of the new parameter 5 given by 5 := h/[l + /i^), (13.870 can be 
rewritten in a slightly more compact form. 



. , . 



1 - 25g„+ig„ + a/ 1 - 4(5g„+ig„ - 45^ (l + g^^^ + ql) 
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Note that if /i G M (or -1 < /i < 1), then -1/2 <6 < 1/2. 

At the "threshold" values of h, h = ±1, corresponding to 6 = ±1/2, we 
can extract the square root in f l3.87p to obtain a rational mapping. This 
is connected with the fact that in the cases h = ±1, the discrete-time sys- 
tem (13.841) under the boundary conditions (13.831) has the trivial solution 
Qn = TQn-i, A„ = 1/(1 + q^-i) (cf. §4.6 in [1]). Thus, the discriminant of 
the quadratic equation (I3.85P vanishes at h = ±1, and the two solutions in- 
tersect. To obtain a nontrivial mapping from this observation, we set = +1 
in (I3.84P and replace g„ with iwn, namely. 



Wn - Wn = (l - W„)A„+i (Wn+l - Wn-l) 

1+1- 

Moreover, we generalize the boundary conditions (I3.83P as 



(1 - W^)An = (1 - wl)An, 



(3.89) 



lim Wn = 7, lim A„ = , } , , 7 -1 

n^±oo n^±oo I 1 -|- 



and assume that \wn — 7I is sufficiently small. Thus, the additional condition 
7 7^ excludes the trivial time evolution u;„ = —Wn-i- Following the same 
procedure as above, we obtain the quadratic equation in A„ that can be 
factorized as 

[(1 + Wn) (1 + Wn-l) An ~ 1] [(1 - ^n) (1 " W„_i) A„ - 1] = 0. 

Substituting its proper solution A„ = 1/ [(1 + Wn) (1 + Wn-i)] into (13.890 . 
we obtain a single equation {wn — l)(u;„+i + 1) = {wn — l){wn-i + 1). By 
a trivial change of the dependent variable Wn ='■ 1 + 2z/?/„ {u ^ 0), we obtain 
the well-known "discrete Volterra equation" [17,49], 

y„(l + l^Vn+l) = 2/n(l + ^Vn-l)- (3.90) 

Therefore, the discrete-time equation (13.901) belongs to the modified Volterra 
hierarchy and not the original Volterra hierarchy [1]; this corresponds to the 
special case where the quadratic equation for the auxiliary variable A„ is 
factorized to provide a rational solution. 

Ultradiscretization [12,13]. We present an intuitively plausible derivation 
of an ultradiscrete analogue of the time-discretized modified Volterra lattice, 
although this may not be a unique ultradiscretization. For the forward time 
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evolution, we rewrite fl3.88p as 



= ^ 

1 + eqn+lQn + 

qn\l + e (g„ + 1/ g„) + V% - equ+i/qn \ 

= r= (3.91) 

1 + eqn+iqn + V 

where e := —25 and F„ := 1 + 2eg„4.ig„ — (l + q^^^ + g^). Note that the 
equation g„ = g„ in the trivial case e = preserves the sign of the depen- 
dent variable on each lattice site n G Z. Thus, we consider a class of solu- 
tions in the limit e — >■ -|-0 such that the positivity of the dependent variable, 
g„ > 0, Vn G Z, is preserved in the nontrivial time evolution. We intro- 
duce the parametrization e =: exp(— L/e), L > and g„_m =: P exp(Q„ m/e), 
p > 0, where m G Z is the discrete time coordinate, which is usually sup- 
pressed. The scaling parameter p can depend weakly on n, m, and £, but 
for brevity, we treat it as a constant. Once this parametrization is substi- 
tuted into (13.911) . the solution {Qn,m} of the initial- value problem depends 
on e via the time evolution. We hypothesize that by choosing p appropri- 
ately, the e-dependence of Q„ in the considered class of solutions becomes 
negligible as e ^ -|-0, and Qn is nonnegative, Qn > 0. Probably, the simplest 
way to justify the latter condition is to modify the zero boundary condi- 
tions (cf. fl3.83p ) to nonzero boundary conditions, lim„^-|-oo = p, and to 
assume g„ > p, Vn G Z. We are not interested in tracing infinitely long tails 
of solitons that decay exponentially as n — ±oo; rather, we prefer to extract 
solitons with compact support in the limit e —>■ +0. This is why we impose 
the zero boundary conditions lim„^-|-oo Qn = on the new variable Qn- 

The condition Yn > 0, guaranteed by the nonnegativity of the discrimi- 
nant of fl3.85p . requires that 

eqn < Qn+i + ^(1 - e2) (l + q^^^) and eg„+i < g„ + ^(1 - e^) (1 + q^) . 

Using the parametrization e = exp(— L/e) and qn = P Gxp{Qn/e) with Qn > 0, 
we can interpret these conditions in the hmit e ^ +0 as —L < Qn+i — Qn < L. 
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By taking the logarithm, (13.911) becomes 



Qn = Qn + e log 



1 + e- 



2Qn 



p +e £ 



% - e- 



Qn + l~Qri—L' 



elog 



Yn = 1 — e £+pe 



2. . / 2 _e -e 



(3.92a) 



(3.92b) 



Taking the aforementioned assumptions into account, fl3.92p reduces in the 
hmit e +0 to the following equation: 



Qn = Qn + max 1^0, Qn + Qn-l - L, 

- max (0, Qn+i + Qn - L 



Qn+l + Qn ~ L 



(3.93) 



As mentioned above, the boundary conditions lim„^-|-oo Qn = are assumed. 
It is desirable to confirm for each solution that the conditions Qn > and 
~ Qn\ < L hold true in the time evolution. Note that the parameter L 
may depend on the discrete time. It is hoped that the issue of the commuta- 
tivity of the ultradiscrete flows defined by (13.931) for distinct values of L will 
be investigated. In this regard, we expect that (I3.93P and its time reversal 
(see below) will comprise an ultradiscrete analogue of the mKdV hierarchy, 
wherein the parameter L labels each flow in the hierarchy. The ultradiscrete 
equation (I3.93P is "linear" in the sense that it is invariant under the rescal- 
ing Qn ^— > kQn, L ^-^ kL, k > 0. Thus, if L is a time-independent constant, 
it is possible to fix L at unity. However, this normalization often changes an 
integer-valued Qn to a fractional value; thus, we rather prefer to leave L as a 
free parameter. Our derivation of the ultradiscrete modified Volterra lattice 
is more or less intuitive and is not mathematically rigorous. 



(I3.93P is more or less intuitive and is not mathematically rigorous. A more 
detailed treatment of all the terms in the numerator of (I3.9ip may lead to a 
more complicated ultradiscrete equation, but we prefer the relatively simple 
equation (I3.93p . Fortunately, for the specific examples that we considered, 
(I3.93P allows the "stable" propagation of solitons and their elastic (but non- 
trivial) collisions; thus, (at least some of) the integrability properties appear 
to be retained in this ultradiscretization. 

The backward time evolution of the time-discretized modified Volterra 
lattice is obtained from fl3.88p through the combined space and time reflection 
n — > —n, m — > —m. Thus, its ultradiscrete analogue is obtained from (I3.93P 
in the same manner. 
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3.5 The lattice Heisenberg ferromagnet model 

The lattice Heisenberg ferromagnet model was proposed in 1982 by several 
different authors [50-53]. The equation of motion for this semi-discrete model 
can be derived almost systematically using either Ishimori's approach [51] 
based on a gauge transformation from the Ablowitz-Ladik lattice or the r- 
matrix formalism based on the fundamental Poisson bracket relation [50,53]. 
However, for a concise and easy-to-understand derivation, we use a more 
heuristic approach based on the zero-curvature representation. 

We start with the Lax pair of the following form in the semi-discrete case: 

Ln = I + XSn, (3.94a) 

Mn = Y^AnLn, (3.94b) 

where the conditions (5'„)^ = / and AnSn = Sn~iAn are assumed. The latter 
condition guarantees the useful relation A„L„ = Ln~iAn. Thus, substituting 
the Lax pair f l3.94p into the zero-curvature condition (12.21) . we obtain 

{Ln^)t + Y^iAn+l - An) = O. 

Noting the identity (/ + AS'„)(/ — XSn) = (1 — A^)/, this results in 

Sn,t = - An. (3.95) 

Because (5'„)^ = / and A„S'„ = Sn-iAn, (13.951) implies the relation 

A-n+l{Sn+l + Sn) = {Sn + >S'„_i)y4„. 

Thus, we choose An as 

An = 2iaSn-l{Sn + Sn-l)'' + 26(^„ + Sn-l)'', 

SO that the above relation is automatically satisfied. Here, a and b are n- 
independent scalars, but they may depend on the time variable t. Substitut- 
ing this expression for An into (I3.95p . we obtain 

Sn,t = A+ [2iaSn-i{Sn + Sn-i)'' + 26(5„ + Sn-i)-'] , (3.96) 

where is the forward difference operator in the discrete space (cf. (12.41) ). 
Note that the evolution equation (13.961) is consistent with the condition 
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{SnY = I for a general I x / matrix Sn- We now consider the simplest non- 
trivial case of Z = 2 and express Sn in terms of the Pauli matrices as 



Sn 



Y,Sli^a,{=:<T.Sn) 

q(3) r.(l) • q(2) 

On iJn '-iJn 

, ; 0(2) r.(3) 
On ~T i-On On 



(3.97) 



Here, Sn = [Sn\ Sn\ Sn^) is a unit vector, i.e., {Sn,Sn) = l- Because 
2{Sn + Sn-i)-' = (Sn + + 5„_i)) [51], dSinD reduces to a sin- 

gle vector equation involving both the scalar product and the vector product. 



Sn + S 



n— 1 



Sn X Sn-1 

1 + {S 

ni Sn—1 71) Sn—l) 



{Sny Sn) 



(3.98) 



The case 6 = gives the lattice Heisenberg ferromagnet model [51,52], while 
the case a = corresponds to its simplest higher symmetry [54]. 

Let us examine the discrete-time case. We mainly consider the time 
discretization of the reduced system (13.981) . and not the general matrix system 



(I3.96p . because the latter problem is expected to be too complicated. We 
start with the Lax matrix Vn of the following form: 



Vn = I + h 



A 

1 _ ^2 " 



(3.99) 



where the conditions (5'„)^ = / and AnSn = Sn-iAn are assumed. The latter 
condition guarantees the useful relation v4„L„ = Ln-iAn- Thus, substituting 
the Lax pair, (I3.94al) and (13.991) . into the zero-curvature condition (12. 6p . we 
obtain 

l{L-' - L-') + - An) = O, 

or equivalently. 



h 



{Sn ~ Sn) — An+l — 



(3.100) 



Under the condition (Sn)"^ = I, the relation AnSn = Sn-iAn is automatically 
satisfied if An takes the general form An = BnSn + Sn^iBn- However, in 
analogy with the semi-discrete case, we employ a more specific form of An, 



An — 2ianSn-l{^Sn + Sn-l) + 2bn[Sn + Sn-l) 



-1 
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which also satisfies the relation AnSn = Sn-iAn. Here, the scalar unknowns 
a„ and 6„ are auxiliary variables. Substituting this expression for An into 
(Kim , we obtain 



h 



( Sn Snj 



2ia„^„_i(S„ + S„_i) ' + 2bn{Sn + Sn^i) ' ■ (S.lOi: 



It only remains necessary to fix the auxiliary variables a„ and 6„. To this end, 
we recall that the time evolution determined by (13.1011) has to be consistent 
with the condition (S'„)^ = /. In the following, we focus on the case of the 
2x2 matrix Sn given by fl3.97p . The requirement tiSn = results in the 
n-independence of a„, thus we set a„ = a. Therefore, fl3.1Ull) reduces to 



^ ( Sn S'n) — 



a h On 

l + (5 ) l + (5„,5„_i) 



(3.102) 



To fix bn, we invoke the procedure presented in section [21 we assume the 
boundary conditions 



lim (5„,5„_i) = 1, lim bn = b. 



(3.103) 



Because {Sn, Sn) = 1, we have det L„ = 1 — (cf. (I3.94ap and (13.971) ). Thus, 
the equality (12.111) derived from (12.101) implies that the determinant of Vn, 



det Vn = det{ I + h 



A 



^ yr)-''-n^n 

det (/ - A^n + hXAn) det 



is an ra-independent quantity. Therefore, both the trace and the determinant 
of Sn — hAn must be ra-independent. The ra-independence of tr(S'„ — hAn) 
is already satisfied by setting a„ = a. The determinant of Sn — hAn can be 
computed as 



det {Sn - hAn) = det Sn - 2ihaSn-i {Sn + Sn-i) ^ - 2/i6„(S'„ + S'„_i) 



det 



1^1 — 2hbn)I + SnSn-i — 2ihaSn-i 



det {Sn + Sn-l) 
(1 + {Sn, Sn^l) - 2hbnY + ^ + ^(H' " (-^n, Sn-l)' 
-2{l + {Sn,Sn-l)) 
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which coincides with its boundary value determined by fl3.103p . This results 
in a quadratic equation in hbn, i.e., 

2{hb^y-2{l + {S^,Sn.l))hbrr 

+ 2{haf + [hh{2 - hh) - [haf] (l + = 0. 

For sufficiently small /i (0 < <^ 1), the proper solution of this quadratic 
equation is given by 



hhn 



2{haf + [hh{2 - hh) - [haf] (1 + g^) 



l + gn + V(l + Qn? - 4(/ia)2 - 2 [hh{2 - hh) - [haf] (1 + gn) 



where gn '■= {Sn, Sn-i)- Substituting this local expression for hhn into fl3.102p . 
we obtain an integrable time discretization of (13.981) : this time discretization 
is essentially equivalent to (B.18) in [15] (see also (26) in [16] and (3.10) 
in [17]). We write the equation of motion for three important cases: the case 
6 = 0, 



ha- 



Sn X Sn-1 
1 + {Sn, Sn-l) 



+ 



1 - 



l-2(/ia)2- 



1 {Sn, Sn-l) 



{l + {Sn,Sn-l)y 

the case hh{2 — hh) = (ha)'^, 

Sn X Sn—l 



Sn + Sn-1 



'n '~-'n 



ha- 



1 + {Sn, Sn^l) 



l-Jl 



and the case a = 0, 
1 



4:{ha) 



{l + {Sn,Sn-l)y 



Sn + Sn-1 



■ ( Sn Sn^ 



Sn-1 



{Sn, Sn) — 1; 



{Sn, Sn) — 1; 



_1 + {Sn, Sn-l) + V (1 + (^n, Sn-l)) " 26 {l + (5„,, 

6:=hh{2-hh), {Sn,Sn) = l. 
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The first and second cases provide time discretizations of tlie lattice Heisen- 
berg ferromagnet model [ fl3.98p with 6 = 0], while the third case gives a 
discrete-time analogue of its simplest higher symmetry [f l3.98p with a = 0]. 
Note that the first case 6 = can be identified with (6.8) in [15]. All these 
difference schemes are seemingly implicit, that is, the value Sn is not written 
explicitly in a closed form in terms of Sn, Sn+i, and Sn-i, but this is not a 
serious drawback. Indeed, we can obtain the power series expansion of Sn in 
h (or 6) to any order successively. Moreover, we can compute the exact value 
of Sn by the following steps. First, we take the scalar product J) et ween the 
equation of motion and Sn+i to obtain an equation for {Sn+i, Sn), wherein 
Sn appears only through the form {Sn+i, Sn)- Then, the proper solution of 
this equation can be found in an 0{h) (or 0{6)) neighborhood of {Sn+i, Sn)- 
Substituting it back into the original equation of motion, we arrive at a linear 
equation for Sn that can be solved straightforwardly. 

In this subsection, we have obtained the time discretizations of the first 
two flows of the lattice Heisenberg ferromagnet hierarchy in a unified way. 
Our derivation is original and easy to follow, but not fully systematic, as 
some parts are based on heuristic treatments. It would be interesting to 
investigate whether these time discretizations can be derived in a systematic 
manner from the time discretizations of the Ablowitz-Ladik lattice obtained 
in subsection 13.21 using Ishimori's approach [51]. 



4 Concluding remarks 

In this paper, we have developed an effective method for obtaining proper 
time discretizations of integrable lattice systems in 1 + 1 dimensions. This 
method, which is based on the zero-curvature condition (12. 6p . allows us to ob- 
tain local equations of motion that can determine the time evolution uniquely 
and explicitly. Using this method, we constructed new time discretizations 
of the Toda lattice, the Ablowitz-Ladik lattice, the Volterra lattice, and 
the modified Volterra lattice, while we obtained the same time discretiza- 
tions of the lattice Heisenberg ferromagnet model and its symmetry as those 
in [15-17]. It should be stressed that this is a systematic method and also 
applies to other integrable lattice systems. As a bonus, we were able to de- 
rive ultradiscrete analogues of the Volterra lattice and the modified Volterra 
lattice involving one arbitrary parameter L, namely, fl3.77p and (13.930 . Each 
of these ultradiscrete equations, as well as its time reversal, appears to form 
a hierarchy of mutually commuting flows labelled by the parameter L. Note 
that the ultradiscrete equations such as (13.770 and (13.930 allow straightfor- 
ward "multicomponent" generalizations. Indeed, for example, if we express 
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the dependent variable and the parameter in the form Un = Un + Un i, 
Ui^^ G Z and L = L^^) + L^^H > 0, L^-') G Z, where i is an irrational num- 
ber, and substitute them into (I3.77p . we can uniquely determine the time 
evolution of the "two- component" system starting from a given initial condi- 
tion. It would be very interesting to investigate the special case wherein the 
irrational i is arbitrarily close to 1. 

A notable feature of our method is that the time discretization does not 
modify the integrable hierarchy to which the original lattice system belongs. 
Thus, the integrals of motion and the functional form of the solutions remain 
invariant; only the time dependence of certain parameters corresponding to 
the angle variables in the solutions is changed (cf. [2,3,18]). Such a time 
discretization can be identified with the spatial part of an auto-Backlund 
transformation of the continuous-time lattice system or, from a more unified 
point of view, an auto-Backlund transformation of the whole hierarchy of 
commuting flows. Therefore, any time discretization obtained by our method 
preserves the qualitative properties of the original integrable hierarchy and 
thus is expected to serve as an excellent scheme for numerical integration 
of the continuous-time lattice system. This is in contrast with other known 
methods that generally modify the integrable hierarchy to find a time dis- 
cretization in the local form. As is illustrated in subsection 13.21 our method 
can also provide "almost integrable" numerical schemes that far surpass the 
Ablowitz-Ladik-Taha local schemes [2,3,5,41] in approximation accuracy. 

In our approach, a local-in-space time discretization is always derived 
from a global-in-space time discretization involving "nonlocal" auxiliary vari- 
ables. The obtained local equations of motion generally determine both the 
forward and backward time evolution uniquely and explicitly under the spec- 
ified boundary conditions; thus, they indeed give a discrete-time analogue of 
evolutionary differential-difference equations. One might consider that the 
original nonlocal time discretization no longer has any use once the local 
discretization has been derived from it; however, this is often not the case. 
In actually solving an initial- value problem, the nonlocal time discretizations 
often provide crucial information on the attributes of the solution, such as 
the real-valuedness, positivity, and rationality with respect to the initial data 
and the parameters; it appears that the local discretizations are not useful 
for establishing such properties directly. Let us illustrate this point with one 
instructive example, namely, the discrete-time Ablowitz-Ladik lattice in the 
case a = b = studied in subsection [321 The global scheme (13.281) (cf. (I3.27P ) 
is free from irrational functions; consequently, under the decaying boundary 
conditions (cf. (13.241) ). the discrete-time updates of the dependent variables 
are given by a rational mapping. On the other hand, the local scheme (I3.34p 
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with fl3.32p involves a square root; thus, it is extremely difficult to foresee 
that the time evolution can be described by a rational mapping. However, 
these two schemes naturally define the same time evolution. Therefore, the 
quantity inside the square root in f l3.32p is always equal to the square of a 
rational function of the dependent variables at the previous moment and the 
parameters; typically, one encounters the form 

where c and d depend on the "step size" parameter h through (13.301) . The 
rational function /„ involves the dependent variables and the parameters 
appearing in (13.341) . and is of order 0{h). In addition, it vanishes when the 
dependent variables are identically zero. For sufficiently small |/?^|(^ 1), this 
square root is extracted as (1 — cd){l + fn)- However, for relatively large \h\, 
this is not obvious; even when cd < 1, the values of the dependent variables 
have to satisfy rather restrictive conditions. One simple way to bypass this 
sign problem is to consider the square root as being defined by its Taylor 
expansion for small \h\ or small amplitudes of the dependent variables and 
then to extend the domain of the definition by analytic continuation; this 
is briefiy explained in subsection 13.21 Thus, this example also illustrates 
the difficulty of the sign problem in computing updates of the dependent 
variables using the local equations of motion. 

In this paper, we mainly considered the "constant" boundary conditions 
at spatial infinity (cf. (12. 9p ) and assumed simple boundary values of the 
dependent variables. Note, however, that our method is not sensitive to 
the boundary values of the dependent variables and is also applicable to 
other boundary conditions, including periodic, nonvanishing, or nonconstant 
boundary conditions. Indeed, as is clear from the description in subsec- 
tion 12.31 one can freely modify the boundary conditions for Vn as long as 
they determine a definite value for the right-hand side of (I2.14p or (I2.15P 
that is ra-independent. Alternatively, one can first set each value of (I2.14p 
or (I2.15P as a "constant" free parameter and then elicit the corresponding 
boundary conditions for the dependent variables. In either case, one should 
take care to identify the proper solution of the resulting algebraic system 
(cf. subsection 12. 4p . We can also understand in a more intuitive way that 
the time discretizations derived in this paper are integrable under other suit- 
able boundary conditions; note that the conservation laws should be derived 
from the local equations of motion using only local operations, i.e., without 
referring to the boundary conditions. 

Very recently, Adler, Bobenko, and Suris [36,37] successfully classified dis- 
crete integrable systems on quad-graphs using the notion of three-dimensional 
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consistency [35]. In particular, under some assumptions, they presented a 
short but complete list of one-field integrable equations defined by polyno- 
mial relations of degree one in each of the four arguments. Their results have 
been attracting a lot of interest from researchers; for example, click on the 
CiteBase link at 

http : //arxiv . org/abs/nlin/0202024l . 
The time discretizations of polynomial lattice systems obtained by our method 
generally contain irrational nonlinearity; the nonlinear terms are determined 
through the proper solution of an "ultralocal" algebraic equation of degree 
higher than one. Thus, such time discretizations essentially lie outside the 
class considered by Adler et al. [36,37], although their three-dimensional con- 
sistency can, in principle, be investigated using the approach outlined in the 
last paragraph of subsection 12. 5[ This probably partially explains why their 
list of one-field integrable equations is short. 

Acknowledgments 

The author thanks Professor J. Matsukidaira, Professor K. Nishinari, Profes- 
sor S. Tsujimoto, Professor K. Toda, Professor T. Takenawa, Dr. F. Sugino, 
Dr. Z. Tsuboi, and Dr. K. Maruno for their useful comments and discussions. 

References 

[1] Y. B. Suris: The Problem of Integrable Discretization: Hamiltonian Ap- 
proach (Birkhauser, Basel, 2003). 

[2] M. J. Ablowitz and J. F. Ladik: Stud. Appl. Math. 55 (1976) 213. 

[3] M. J. Ablowitz and J. F. Ladik: Stud. Appl. Math. 57 (1977) 1. 

[4] M. J. Ablowitz and J. F. Ladik: J. Math. Phys. 17 (1976) 1011. 

[5] T. R. Taha and M. J. Ablowitz: J. Comp. Phys. 55 (1984) 192. 

[6] Y. B. Suris: Phys. Lett. A 234 (1997) 91. 

[7] Y. B. Suris: Inv. Probl. 13 (1997) 1121. 

[8] R. Hirota: The Direct Method in Soliton Theory (Cambridge Univ. 
Press, Cambridge, 2004) edited and translated by A. Nagai, J. Nimmo 
and C. Gilson. 



58 



[9] H. Flaschka: Phys. Rev. B 9 (1974) 1924. 

[10] H. Flaschka: Prog. Theor. Phys. 51 (1974) 703. 

[11] S. V. Manakov: Sov. Phys. JETP 40 (1975) 269. 

[12] D. Takahashi and J. Matsukidaira: Phys. Lett. A 209 (1995) 184. 

[13] T. Tokihiro, D. Takahashi, J. Matsukidaira and J. Satsuma: Phys. Rev. 
Lett. 76 (1996) 3247. 

[14] F. W. Nijhoff, G. R. W. Quispel and H. W. Capel: Phys. Lett. A 95 
(1983) 273. 

[15] G. R. W. Quispel, F. W. Nijhoff, H. W. Capel and J. van der Linden: 
Physica A 125 (1984) 344. 

[16] F. W. Nijhoff, H. W. Capel, G. L. Wiersma and G. R. W. Quispel: Phys. 
Lett. A 103 (1984) 293. 

[17] F. Nijhoff and H. Capel: Acta Appl. Math. 39 (1995) 133. 

[18] F. Kako and N. Mugibayashi: Prog. Theor. Phys. 61 (1979) 776. 

[19] Y. B. Suris: Inv. Probl. 16 (2000) 1071. 

[20] S. J. Orfanidis: Phys. Rev. D 18 (1978) 3822. 

[21] M. J. Ablowitz: Stud. Appl. Math. 58 (1978) 17. 

[22] K. Kajiwara and A. Mukaihira: J. Phys. A: Math. Gen. 38 (2005) 6363. 
[23] S. J. Orfanidis: Phys. Rev. D 18 (1978) 3828. 

[24] E. Date, M. Jimbo and T. Miwa: J. Phys. Soc. Jpn. 51 (1982) 4116. 
[25] T. Sadakane: J. Phys. A: Math. Gen. 36 (2003) 87. 
[26] G. L. Lamb, Jr.: Rev. Mod. Phys. 43 (1971) 99. 

[27] A. C. Scott, F. Y. F. Chu and D. W. McLaughlin: Proc. IEEE 61 (1973) 
1443. 

[28] R. M. Miura (editor): Bdcklund Transformations, the Inverse Scattering 
Method, Solitons, and Their Applications (Lect. Notes in Math. 515, 
Springer, Berlin, 1976). 



59 



[29] B. G. Konopelchenko: Phys. Lett. A 87 (1982) 445. 

[30] H. D. Wahlquist and F. B. Estabrook: Phys. Rev. Lett. 31 (1973) 1386. 

[31] V. Papageorgiou, B. Grammaticos and A. Ramani: Phys. Lett. A 179 
(1993) 111. 

[32] V. P. Spiridonov, S. Tsujimoto and A. S. Zhedanov: Gommun. Math. 
Phys. 272 (2007) 139. 

[33] K. Kajiwara and Y. Ohta: J. Phys. Soc. Jpn. 77 (2008) 054004. 

[34] R. Hirota: J. Phys. Soc. Jpn. 43 (1977) 2079. 

[35] F. W. Nijhoff and A. J. Walker: Glasgow Math. J. 43A (2001) 109. 

[36] V. E. Adler, A. I. Bobenko and Y. B. Suris: Gommun. Math. Phys. 233 
(2003) 513. 

[37] V. E. Adler, A. I. Bobenko and Y. B. Suris: Funct. Anal. Appl. 43 
(2009) 3. 

[38] R. K. Dodd: J. Phys. A: Math. Gen. 11 (1978) 81. 
[39] Y. B. Suris: Phys. Lett. A 206 (1995) 153. 

[40] F. Pempinelli, M. Boiti and J. Leon: Proceedings of the First Workshop 
on Nonlinear Physics, Theory and Experiment (Gallipoli, Italy, 1995) 
edited by E. Alfinito et al. (World Scientific, Singapore, 1996) pp. 261. 

[41] T. R. Taha and M. J. Ablowitz: J. Gomp. Phys. 55 (1984) 203. 

[42] L. D. Faddeev and L. A. Takhtajan: Hamiltonian Methods in the Theory 
of Solitons (Springer, Berlin, 1987) translated by A. G. Reyman. 

[43] J. Matsukidaira, J. Satsuma, D. Takahashi, T. Tokihiro and M. Torii: 
Phys. Lett. A 225 (1997) 287. 

[44] D. Takahashi and J. Matsukidaira: J. Phys. A: Math. Gen. 30 (1997) 
L733. 

[45] D. Takahashi and J. Satsuma: J. Phys. Soc. Jpn. 59 (1990) 3514. 
[46] R. Hirota: J. Phys. Soc. Jpn. 35 (1973) 289. 
[47] M. Wadati: Prog. Theor. Phys. 59 (1976) 36. 



60 



[48] M. J. Ablowitz and J. F. Ladik: J. Math. Phys. 16 (1975) 598. 

[49] R. Hirota and S. Tsujimoto: J. Phys. Soc. Jpn. 64 (1995) 3125. 

[50] E. K. Sklyanin: Punct. Anal AppL 16 (1982) 263. 

[51] Y. Ishimori: J. Phys. Soc. Jpn. 51 (1982) 3417. 

[52] F. D. M. Haldane: J. Phys. C: Sohd State Phys. 15 (1982) L1309. 

[53] L. Faddeev: Recent Advances in Field Theory and Statistical Mechanics 
(Lcs Houchcs Summer School, Session XXXIX, 1982) edited by J.-B. 
Zuber and R. Stora (North-Holland, Amsterdam, 1984) pp. 561. 

[54] V. S. Gerdjikov, M. I. Ivanov and Y. S. Vaklev: Inv. Probl. 2 (1986) 
413. 



61 



